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CHAPTER  1 
INTRODUCTION 


In  this  study  some  numerical  issues  related  to  the  computational 
solution  of  a  generalized  form  of  the  algebraic  matrix  Rlccatl  equation 
are  examined.  The  approach  herein  utilizes  the  generalized  eigenprob- 
lem  formulation,  which  provides  a  powerful  framework  for  the  solution 
of  quite  general  forms  of  algebraic  Rlccatl  equations  arising  in  both 
continuous-  and  discrete-time  applications.  This  general  form  is 
derived  from  control  and  filtering  problems  for  systems  in  generalized 
(or  implicit  or  descriptor)  state  space  form.  These  equations  play 
fundamental  roles  in  the  analysis,  synthesis,  and  design  of  linear- 
quadrat  lc-Gausslan  control  and  estimation  systems  as  well  as  in  other 
areas  of  applied  mathematics.  A  representative  sample  of  applications 
may  be  found  in  [  1  ]— [  A] .  It  is  not  our  purpose  here  to  survey  the 
extensive  literature  available  on  Rlccatl  equations,  but,  rather  we 
refer  the  reader  to,  for  example  [l]-[4]  for  references. 

The  method  exploited  here  is  a  variant  of  the  classical  eigen¬ 
vector  approach  to  Rlccatl  equations.  Martensson  [5]  is  one  of  the 
beat  suunaries  of  the  eigenvector  approach  to  solving  algebraic  Rlccatl 
equations.  However,  the  use  of  eigenvectors  often  encounters  numerical 
difficulties  in  practical  computation,  especially  when  the  correspond 
ing  eigenvalues  are  closely  spaced.  Thus,  the  method  to  be  preferred 
and  the  basis  of  this  work  employ  Schur  vectors  instead  of  eigenvec¬ 
tors  because  Schur  vectors  are  more  reliably  and  accurately  computed. 
The  Schur  method  was  first  examined  in  detail  by  Laub  [ 6  ] ,  [  7  ]  and 
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then  extended  by  Pappas,  et.al.  [8]  and  Emami-Naeini  [9]  for  the 
discrete-time  problem,  and  by  Van  Dooren  [10]  for  the  continuous-  and 
discrete-time  problem.  The  very  general  formulation  herein  was  derived 
by  Laub  [11]  and  expanded  by  Lee  [12]. 

The  numerical  Issue  first  examined  In  this  thesis  Is  the  deriva¬ 
tion  of  a  Newton-type  iterative-refinement  procedure  for  the  general¬ 
ized  problem  formulation.  Klelnman  [ 1 3 ]  and  Hewer  [14]  derived  methods 
of  this  type  for  the  standard  Rlccatl  equation  In  the  continuous-  and 
discrete-time  characterizations,  respectively. 

Perhaps  the  most  important  issue  studied  here  is  that  of  condi¬ 
tioning  of  the  Rlccatl  problem,  that  Is,  to  define  a  measure  of  the 
sensitivity  of  the  numerical  solution  to  changes  in  the  data.  Although 
results  are  available  on  conditioning  for  problems  such  as  linear 
equations  [ 15 ] ,  [  16]  and  eigenvalues  [l7]-[l9],  little  is  published  on 
the  conditioning  of  the  Rlccatl  problem.  A  recent  work  In  this  area 
for  the  continuous-time  problem  is  by  Byers  [  20] .  The  results  of 
numerical  experiments  examining  various  condition  measures  are  reported 
herein. 

Since  the  behavior  of  many  physical  systems  in  engineering  Is 
first  modeled  by  second-order  differential  equations  with  very  special 
properties  (i.e.,  symmetry  and  definiteness),  we  examine  exploiting 
this  structure  in  certain  computations  of  interest.  Namely,  we  exploit 
the  structure  in  solving  associated  Rlccatl  problems  and  In  testing  for 
controllability  and  observability  of  the  second-order  models. 

4 


1.1  Main  Contributions  of  the  Thesis 


We  regard  the  derivation  of  condition  measures  for  the  continuous- 
and  discrete-time  generalized  algebraic  Rlccati  equation  and  the 
evaluation  of  their  numerical  behavior  in  special  situations  as  the 
main  contributions  of  this  thesis.  More  specifically,  we  may  list  the 
following  contributions. 

1.  The  derivation  of  a  Newton-type  iterative-refinement  procedure 
for  the  continuous  and  discrete  algebraic  Rlccati  equation  that 
has  monotonic  convergence  that  is  quadratic  in  the  vicinity  of 
the  solution. 

2.  Establishing  the  equivalence  between  the  given  generalized 
optimal  control  problem  solution  and  the  solution  of  a  more 
standard  optimal  control  problem  in  certain  cases. 

3.  Derivation  of  condition  numbers  for  the  solution  of  the  gener¬ 
alized  Rlccati  equation. 

4.  A  new  algorithm  for  applying  a  change  of  coordinate  transforma¬ 
tion  to  the  system  model  when  the  Rlccati  problem  is  solved 
using  the  generalized  eigenproblem  approach. 

5.  Numerical  experimentation  to  illustrate  the  behavior  of  known 
and  speculated  condition  estimates  for  the  Rlccati  problem. 
The  numerical  results  demonstrate  that  a  single  satisfactory 
condition  estimate  is  not  available. 

6.  Exploitation  of  the  structure  of  second -order  models  to  solve 
efficiently  the  velocity  feedback  optimal  control  problem  using 
a  generalized  Rlccati  equation. 
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7.  Formulation  of  useful  tests  for  controllability  and  observabil¬ 
ity  of  second -order  models  directly  in  terms  of  conditions  on 
the  model  matrices. 

8.  Development  of  a  portable  FORTRAN  software  package  for  the 
efficient,  reliable  solution  of  generalized  algebraic  Riccatl 
equations. 

1.2  Outline  of  Chapters 

In  Chapter  2  we  review  some  important  concepts  of  numerical 
analysis,  which  are  used  heavily  in  the  succeeding  analysis.  The 
concepts  of  numerical  stability  and  conditioning  are  presented,  and 
there  is  a  brief  discussion  of  first-order  perturbation  analysis.  The 
application  of  Schur  techniques  for  the  solution  of  generalized  alge¬ 
braic  Riccatl  equations  is  reviewed.  Very  general  optimal  control  and 
filtering  problems  are  formulated,  which  result  in  generalized  algebraic 
Riccatl  equations  for  the  continuous-  and  discrete-time  cases. 

Chapter  3  derives  a  Newton-type  iterative  procedure,  which  we 
employ  for  Improving  the  numerical  accuracy  of  the  Schur  solution  of 
generalized  Riccatl  equations  or  for  calculating  new  solutions  when 
problem  parameters  are  changed  by  a  small  amount.  Newton's  method  for 
the  standard  algebraic  Riccatl  equation  in  continuous-  and  discrete¬ 
time  formulations  is  reviewed  first.  Equivalence  between  the  general¬ 
ized  solution  and  a  certain  standard  solution  is  established.  The 
Newton  procedure  is  then  extended  to  the  generalized  case. 

The  subject  of  conditioning  of  the  generalized  algebraic  Riccatl 
equation  is  examiner  C  .pter  4.  The  desired  form  of  condition 
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estimate  Is  defined.  Previous  work  In  the  area  of  conditioning  of  the 
Rlccstl  problem  is  reviewed.  New  condition  estimates  are  then  derived 
using  a  first -order  perturbation  analysis.  Balancing  to  improve  the 
condition  of  the  problem  is  next  considered.  Balancing  of  the  general¬ 
ized  elgenproblem  is  discussed  first.  A  new  algorithm  for  applying  a 
change  of  coordinates  transformation  to  the  system  model  when  the 
Riccati  problem  is  solved  in  the  generalized  elgenproblem  formulation 
is  presented. 

Numerical  results  for  the  solution  of  the  generalized  Riccati 
equation  for  special  cases  are  presented  in  Chapter  5.  The  software 
package  developed  as  a  research  tool  to  aid  in  the  studies  of  this 
thesis  is  discussed.  Highlights  of  the  package  capabilities  are 
given.  The  numerical  algorithms  employed  and  sources  for  the  software 
are  discussed  and  referenced.  The  results  of  three  specially  designed 
examples  are  examined  to  illustrate  the  behavior  of  the  Riccati  solu¬ 
tion  and  the  ability  of  the  known  and  new  condition  estimates  to 
predict  the  numerical  accuracy  of  the  solution. 

In  Chapter  6  second -order  models  are  considered.  The  second -order 
model  structure  is  defined  in  the  large  space  structure  framework.  The 
generalized  algebraic  Riccati  equation  is  considered  for  this  problem, 
and  ways  are  explored  to  take  computational  advantage  of  the  second- 
order  formulation  to  solve  the  Riccati  equation.  The  velocity  feedback 
problem  is  shown  to  reduce  to  a  simple  form  in  this  framework.  Tests 
for  controllability  (stabilizabillty)  and  observability  (detectability) 


are  then  derived  in  terms  of  the  model  matrices  for  second  order 
models. 

Finally  in  Chapter  7  we  state  the  conclusions  we  have  drawn  from 
this  work  and  make  recommendations  for  future  research  in  this  area. 

An  appendix  is  included  which  contains  the  description  (from  the 
actual  software)  of  all  the  FORTRAN  subroutines  accumulated  in  the 
software  package. 
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CHAPTER  2 


BACKGROUND 


la  this  chapter  we  review  some  Important  concepts  and  results 
which  form  the  basis  for  the  work  In  the  succeeding  chapters.  The 
first  section  discusses  the  numerical  analysis  concepts  of  stability 
and  conditioning.  The  second  section  briefly  presents  the  perturbation 
theory  necessary  In  later  analysis.  The  third  section  reviews  results 
employing  Schur  techniques  for  the  solution  of  algebraic  matrix  Rlccati 
equations.  Very  general  optimal  control  and  filtering  problems  are 
formulated  which  result  in  generalized  algebraic  Rlccati  equations 
(GARE)  for  the  continuous-  and  discrete-time  cases.  Solutions  for 
these  GARE  are  stated  In  terms  of  a  generalized  eigenproblem. 

2.1  Numerical  Stability  and  Conditioning 

We  shall  now  Introduce  the  important  concepts  of  numerical 
stability  and  conditioning.  The  following  framework  will  be  useful  for 
the  understanding  of  these  concepts: 


(2.1) 


problem 

or 

model 


solutions 


Give  d  e  D,  we  want  to  compute  f(d)  e  S.  However,  frequently  only  d* 
(near  d)  is  known  and  only  f*  (an  algorithm  to  approximate  f)  is 
available.  Figure  2.1  is  a  geometrical  view  of  this  problem-solving 
process.  Therefore,  in  this  framework  we  seek  f(d),  but  compute 
f*(d*). 
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FIGURE  2.1  Geometric  view  of  the  effect  of  data  uncertainty 
and  computational  errors  In  solving  a  problem. 

2.1.1  Numerical  Stability 

Definition  2.1:  An  algorithm  f*  is  numerically  stable  If  for  all  d 
contained  In  D,  there  exists  a  d*  contained  In  D  and  near  d  such  that 
f(d*>  (the  exact  solution  of  a  slightly  perturbed  problem)  is  near 
f*(d)  (the  computed  solution). 

That  Is,  we  expect  that  a  stable  algorithm  will  not  introduce  any 
inaccuracies  Into  the  solution  larger  than  those  present  in  the  data. 

2.1.2  Problem  Conditioning 

Definition  2.2:  If  f(d*)  (the  exact  solution  of  a  slightly  perturbed 
problem)  Is  near  f(d)  (the  true  solution),  the  problem  Is  said  to  be 
well-conditioned.  If  f(d*)  may  potentially  differ  greatly  from  f(d), 
the  problem  is  ill-conditioned. 
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The  four  possible  combinations  of  stability  and  conditioning  are 
shown  in  Figure  2.2. 


a)  Well-conditioned  problem 
Stable  computation 


b)  Well-conditioned  problem 
Unstable  computation 


c)  Ill-conditioned  problem 
Stable  computation 


d)  Ill-conditioned  problem 
Unstable  computation 


FIGURE  2.2  Possible  combinations  of  problem  conditioning  and  computa¬ 
tion  stability  in  the  numerical  solution  of  a  problem. 
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Although  It  nay  be  impossible  for  the  numerical  solution  process 
to  guarantee  an  accurate  answer  to  an  ill-conditioned  problem,  it  is 
desirable  for  the  program  to  recognize  ill -conditioning  of  the  computa¬ 
tions  and  report  this  fact  to  the  user.  Therefore,  it  is  desirable  to 
associate  data  with  a  computing  problem  which  reflects  the  overall 
sensitivity  of  the  solution  to  changes  in  the  data  (l.e.,  a  condition 
number).  He  can  do  this  as  follows  [21]: 

Definition  2.3:  Let  D  and  S  be  finite  dimensional  metric  spaces  with 
metrics  pd  and  p8,  respectively.  Let  f(d)  be  the  computing  problem 
under  consideration.  As  before 
f  :  D  - »  S 

The  absolute  asymptotic  condition  number  is 

p8(f(d),f(d*)) 

ie(f(d))  :■  Him  sup  -  (2.2) 

6-0  pd(d,d*)-«  6 

Relative  condition  can  be  defined  similarly  [21].  Figure  2.3 
Illustrates  the  concept  of  a  condition  number.  One  can  see  that  k  is  a 
measure  of  the  sensitivity  of  the  solution  to  perturbations  in  the 
data,  k  can  be  interpreted  as  the  factor  by  which  the  uncertainty,  6, 
is  multiplied  by  in  the  problem-solving  process.  Condition  numbers  are 
generally  estimated  using  perturbation  theory. 
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FIGURE  2.3  Geometric  interpretation  of  the  condition  number,  K. 


2.1.3  Condition  Number  of  a  Matrix 

Definition  2.4:  Let  1*1  be  a  consistent  matrix  norm  on  Rnxn  satisfy- 


III  -  1 

with  a  consistent  vector  norm  on  Rn.  Let  A  e  Rnxn ,  then 
k(A)  -  IAI  IA~ll  . 


(2.3) 


(2.4) 


Now  consider  the  computation  of  the  Inverse  of  a  slightly  perturbed 
matrix,  i.e.,  find  (A  +  E)_1 . 

Theorem  2.1:  If  IA”1!  •  I  El  <  l,  then 


IA-Ma+E)-1!  < 

I  A-1 1 


< (A)  I  El 
IAI 

1-k(A) IEI 
IAI 


(2.5) 


F*  T*  1 
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The  left  hand  side  of  the  inequality  is  the  relative  error  in 


(A  +  E)_l .  If  E  is  sufficiently  snail,  the  right  side  is  effectively 
k(A)  lEI/lAl .  Therefore,  the  theorem  states  that  the  relative  error  in 
A  +  E  nay  be  magnified  by  as  aich  as  k(A)  in  calculating  (A  +  E)“l. 
For  this  reason,  ic(A)  is  called  "the  condition  number  of  A  with  respect 
to  inversion."  Note  that 

1  -  111  -  IAA“ll  <  IAI  IA'll  -  ic(A)  (2.6) 
2.1.4  Error  Analysis 

In  the  numerical  solution  of  real  problems,  we  compute  f*(d*) 
Instead  of  f(d),  so  we  would  like  to  estimate  lf(d)  -  f*(d*)l.  This  is 
the  basic  goal  of  error  analysis.  There  are  two  main  types  of  error 
analysis  of  numerical  computation  and  are  referred  to  as  forward  and 
backward  error  analysis. 

Ir.  the  forward  error  analysis,  one  attempts  to  obtain  a  bound  on 
the  error  in  the  final  result  by  starting  with  the  original  problem  and 
following,  step  by  step,  the  effect  of  computational  errors  (round-off) 
and  original  data  uncertainty.  However,  the  resulting  bounds  are 
usually  hopelessly  pessimistic, and  the  analysis  itself  is  very  compli¬ 
cated  for  all  but  simple  problems  [16],  [22]. 

In  backward  error  analysis,  one  does  not  attempt  to  compute  lf(d) 
-  f*(d*)(,  but  rather  one  attempts  to  determine  how  close  the  problem 
actually  solved  is  to  the  original  problem.  This  is  illustrated 
graphically  in  Figure  2.4.  In  this  technique,  one  uses  error  bounds  to 
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show  that  the  computed  solution  of  a  given  problem  Is  near  the  exact 
solution  of  a  slightly  perturbed  problem.  This  is  sufficient  to  ensure 
that  the  algorithm  that  did  the  computations  is  stable. 


* 


2.1.5  Role  of  Orthogonal  Matrices 

The  class  of  orthogonal  transformations  has  a  special  role  in 
numerical  computations. 

Definition  2.5:  An  orthogonal  matrix  0  is  a  square  matrix  with  ortho¬ 
normal  columns.  That  is 

U  e  Rnxn,  U  -  (ult  u2,  ...  ,  Un) 

and 


Properties  of  orthogonal  matrices  important  to  numerical  computations 
are 

Theorem  2.2:  Let  U  be  orthogonal,  then 
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I'L-tr""  jr 


l*J 


V 


(2.7) 
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1)  UTU  -  UDT  -  I  ; 

2)  IAUI2  -  IUAI2  -  IAl2. 

Proof :  1)  0TU  -  I  follows  from  the  orthonormality  of  the  columns  of 

U.  It  follows  then  that  U?  Is  the  Inverse  of  U.  Since  by  definition 
a  matrix  commutes  with  its  inverse,  we  have  UU^  -  I.  Note  that  this 
implies  that  the  rows  of  U  are  also 'orthonormal. 

2)  To  prove  this  part,  we  make  use  of  the  facts  that  for  A  e 
R™™,  IAI22  -  IAtAI2  and  IATI2  -  IAl2  ([15],  Theorem  4.2.10). 
Now, 

IUAI22  -  I(UA)TUA12  -  IATUT0AI2  -  IATAI2  -  I Al 22 
and,  therefore, 

IAUI2  -  l(A0)Tl2  -  ,oV,2  -  iati2  -  iai2. 

From  this  it  is  easy  to  see  that  orthogonal  matrices  have  three  advan¬ 
tages.  First,  they  are  easy  to  Invert;  U“l  ■  UT.  Second,  orthogonal 
matrices  are  perfectly  conditioned  with  respect  to  l*l2: 

k(U)  -  IUI2IUtI2  -1.1-1.  (2.8) 

A  third  advantage  of  orthogonal  transformations  is  that  perturbations 
in  the  result  can  be  accounted  for  by  a  perturbation  of  the  same  size 
in  the  original  problem.  For  example,  having  computed  U^AII,  we 
Introduce  an  error  F  into  the  result.  If  we  set  E  -  UFU^,  then  IEI2 
-  t  F 1 2  and 

Ut(A  +  E)U  -  UtAU  +  F  .  (2.9) 
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This  makes  orthogonal  transformations  ideal  for  backward  error  analy¬ 
sis.  That  is,  Multiplications  by  orthogonal  transforaatlons  are 
backward  stable. 

2.2  Perturbation  Analysis 

For  the  purposes  of  this  thesis,  we  eaploy  a  technique  coaaonly 
known  as  first-order-perturbation  theory.  This  technique  reveals  the 
effects  of  perturbations  in  the  problem  on  the  solution.  The  techrlque 
is  generally  applied  in  a  three-step  procedure.  First,  a  fora  is 
chosen  for  the  approximate  problem  solution.  Then  the  approxiaatlon  is 
substituted  into  an  equation,  and  all  terns  second-order  or  higher  in 
small  quantities  are  deleted.  Finally,  the  resulting  linear  equation 
is  solved  for  the  unknown  in  the  approximation. 

To  illustrate  the  technique,  consider  the  problem  of  estimating  (A 
+  E)-1,  where  A  is  nonsingular  and  E  is  a  "small-perturbation"  matrix. 
First  we  chose  a  fora  for  the  approximation  as 

(A  +  E)-1  *  A”l(I  -  H),  (2.10) 
where  H  is  a  "small  perturbation”  matrix.  Substituting  the  approxima¬ 
tion  into  the  equation 

(A  +  E)(A  +  E)_I  -  I,  (2.11) 

we  obtain 

(A  +  E)A“l  (I  -  H)  -  I  -  H  +  EA_l  -  EA“lH  *  I  .  (2.12) 
Dropping  the  term  EA_1H,  which  is  assumed  small  in  comparison  to  E  and 


H,  and  solving  the  resulting  equation  gives  us 
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(A  +  E)_1  «  A~l  (I  -  EA"1).  (2.14) 

We  can  now  loosely  derive  Che  condition  nunber  of  a  matrix  with 
respect  to  Inversion.  From  (2.14)  we  have 

I (A  +  E)"1  -  A-1 I  *  IA-1EA-1I  <  IA"1I2IEI  (2.15) 

Bence,  we  have 


<  .A-‘|.E1  -  .(A)  ill 


IA-‘| 


(2.16) 


IAI 


This  differs  from  the  bound  of  Theorem  2.1  (equation  2.5)  by  the  factor 

1  (2.17) 

1-k(A)|E| 

IAI 

which  la  negligible  when  I  El  Is  small. 

2.3  Solution  of  the  Generalized  Algebraic  Rlccati  Equation  (GARE)  as  a 

Generalized  Eigenproblem 


In  this  section  we  shall  present  the  method  of  solving  the  GARE 
(both  continuous-  and  discrete-time  cases)  via  a  generalized  eigenprob¬ 
lem.  First,  the  GARE  resulting  from  the  optimal  regulator  problem  will 
be  derived.  Then  the  optimal  filtering  problem  will  also  be  stated,  and 
the  corresponding  GARE  given.  Finally,  the  solutions  to  the  GARE  will 
be  formulated  utilizing  the  generalized  real  Schur  form  of  the  general¬ 
ized  eigenproblem. 

2.3.1  Optimal  Regulator  -  Continuous-Time  Problem 

Consider  the  following  general  time-invariant  deterministic  linear 
optimal  regulator  problem: 


System: 


Ex(t)  «  Ax(t)  +  Bu(t)  ;  x(t0)  -  xq 


(2.18) 
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% 

>5 

M 


$ 


<st 


Criterion:  J  -  X  }m  (yTQy  +  uTBu  +  2  xTSu)dt 
^  n 


where 


x  e  Rn;  u  e  R®;  y  e  Rr; 

E,A  e  R“m;  B  e  R11*®;  C  e  Rrxn, 

Also,  0  -  0T  e  R**1;  R  -  RT  e  R®*®;  S  e  Rnx®; 


(2.19) 


(2.20) 


R  >  0  ;  E  nonsingular. 


Application  of  Hamilton -Jacobi  theory  gives  rise  to  the  following 


equations : 


h  *  ^  (xTCQCx  +  2xTSu  +•  uTRu)  +  pT(Ax  +  Bu  -  Ex) 
3h  drdhi  _  .  .  .  _• 

yp“3t^^*°"A*4'Bu*'Ex 


"  at 


{it  }  .  o  .  CTQCx  ♦  Su  +  ATp  +  ETp 


3h  d  r 3h  i  T  _  T 

-0-Sx  +  Ru  +  Bp 


(2.22) 


(2.23) 


(2.24) 


(2.25) 


where  h  Is  the  scalar  Haalltonlan,  and  p(t)  Is  the  costate  vector, 
p  e  Rn. 

In  the  usual  treataent  E  »  I  and  S  ■  0.  equation  (2.25)  is  solved 
for  u,  and  u  la  then  substituted  In  (2.23)  and  (2.24)  to  obtain  the 
Hamiltonian  system 


x  A  -BR“l  1 

•  “  T  T 

pj  [_-C  QC  -A 


(2.26) 


If  one  makes  the  "Rlccatl  substitution’*  p  -  Xx ,  we  are  led  to  the 
algebraic  Rlccatl  equation  for  X 


NWC  TP  6521 


atx  +  xa  -  xbr_1btx  +  ctqc  -  o 


(2.27) 


If  Che  pair  (A,B)  is  stabllizable  and  Che  pair  (C,A)  is  detectable, 
Chen  a  unique  non-negaCive  def iniCe  solucion  X  -  X^  >  0  exists  such 


ChaC  the  linear  control  law 
u°  -  -R-1BTx 


(2.28) 


stabilizes  the  system  and  is  optimalT  in  the  sense  Chat  it  minimizes  the 
criterion  (2.19)  over  all  other  control  laws  [l]. 


In  the  more  general  setting  we  have  the  following  Hamiltonian 


system 


A-BR"lST 


-br-1bt 


sr_1st-ctqc  -(a-br-1st)t 


(2.29) 


If  one  makes  the  "Rlccatl  substitution”  p  *  XEx,  we  are  led  to  the  GARE 

o  -  (a-br-ist)txe  +  etx(a-br-1st)  -  etxbr-ibtxe  +  ctqc  -  sr-1st 


-  ATXE  +  ETXA  -  (BTXE  +  ST)TR-1  (BTXE  +  ST)  +  CTQC  (2.30) 

and  the  optimal  control  law  is  given  by 


HPJ 


u°  -  -R_1(ST  +  BTXE)x 


(2.31) 


To  solve  for  X  «  XT  0,  one  can  use  the  technique  suggested  by 
Pappas,  et.al.  [8]  for  the  discrete  time  ARE  as  expanded  by  Lee  [12]. 
Before  presenting  this  technique,  we  will  formulate  the  other  problems 


of  Interest. 


2.3.2  Optimal  Filter  -  Continuous -Time  Problem 


Consider  the  following  time -Invariant  linear  optimal  observer 


& 


problem 


System:  Ex(t)  -  Ax(t)  +  Bu(t)  +  Dv(t), 


(2.32) 


where 


[«!] 


y(t)  -  Cx(t)  +  w(t ) ,  t  >  t0, 


,Jt+r 


e  R  is  a  white  noise  process  with  intensity 


Jt 

r 


(2.33) 


Also,  (2.20)  applies,  E  is  nonsingular  and  D  e  Rnx&.  Furthermore, 
the  initial  state  x(tg)  is  uncorrelated  with  v  and  w,  and  u(t)  is  a 
given  input  to  the  system. 

The  optimal  observer  is  given  by 

Ei(t)  -  Ax(t)  +  Bu(t )  +  K°[y(t)  -  Cx(t)]  , 
which  minimizes  the  weighted  mean  square  reconstruction  error 

B  {eT(t)We(t)f  , 

where 


(2.34) 


(2.35) 


e(t)  -  x(t)  -  x(t);  W  -  >  0  .  (2.36) 

Let 

S  -  ES  ,  (2.37) 

Then  K°  is  given  by 

K°  -  [S  +  EXCT]  R-1  ,  (2.38) 

where  X  solves  the  following  GARE: 

0  -  (A-SR-lC)XET  +  EX(A-SR-lC)T  -  exctr_1cxet  +  dqdt-  SR_1ST.  (2. 


One  can  see  that  this  GARE  is  the  dual  to  (2.30),  and,  therefore , 
we  can  solve  for  X  by  the  same  technique. 
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Consider  the  following  general  discrete-tine,  tine-invariant 
deterministic  optimal  regulator  problem: 


System:  E*k+1  *  Axk  +  Buk  5  x(0)  *  xq 


(2.40) 


yk  "  Cxk 


Criterion:  J  -  i  f  ^k^k  +  \TRuk  +  2xfcTsuk) 


(2.41) 


where  (2.20)  and  (2.21)  apply  also  except  we  only  require  R  0. 
Application  of  the  discrete  maximum  principle  [23]  gives  rise  to  the 
set  of  equations 

hfc  “  i  xkTCT0Cxk+  uJru^  2xkTSuk  +  pJL  Ax^  Bu^  (2.42) 


E\  -  ^  ■  °T«\  *  s\  * 

3xk 


(2.43) 


T  T 

__  -  0  -  R^  +  S  xk  +  B  pk+1 
3uk 


(2.44) 


3Pk+l 


Axk  +  Buk 


(2.45) 


where  h^  Is  the  scalar  Hamiltonian,  and  p^  e  Rn  is  the  costate 


vector . 


Assuming  for  the  moment  that  R  is  invertible,  we  can  solve  (2.44) 
for  ufc  and  substitute  into  (2.43)  and  (2.45).  Also,  we  make  the 
"Riccati  substitution"  p^  -  XExk  to  obtain  the  following  discrete 
version  of  the  GAR£  (for  nonsingular  X): 
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ETXE  -  (A  -  BR-1ST)T(X-1  +  BR_IBT)  \a  -  BR-lST) 

+  CTQC  -  SR*lST  .  (2.46) 

This  equation  can  be  algebraically  manipulated  Into  the  following 
equivalent  forma: 

T  T 

ETXE  -  (A  -  BR"lST)  X(A  -  SR“lST)  -  (A  -  BR-lST)  XB(R 

+  BTXB)  lBTX(A  -  BR" 1 )  +  CTQC  -  SR_1ST  (2.47) 

-  ATXA  -  (ATXB  +  S)(R  +  BTXB)_l(ATXB  +  S)T  +  CTQC.  (2.48) 

Note  that  (2.47)  does  not  require  X-1  explicitly,  and  (2.48)  does  not 
require  X-1  or  R-1  explicitly. 

Once  X  la  determined  from  the  above  GARS,  the  linear  optimal 
feedback  la  given  by 

u®  -  -(R  +  BTXB)-1  (ATXB  +  S)\  .  (2.49) 

Aa  In  the  continuous -time  case,  stabillzabllity  and  detectability 
assumptions  assure  the  existence  and  uniqueness  of  X  •  X^  0  such 
that  (2.49)  Is  a  stabilizing  feedback. 

2.3.4  Optimal  Filter  -  Discrete-Time  Problem 

Consider  the  following  discrete-time,  time-invariant  linear 
optimal  observer  problem 

System:  Ex^+i  ■  Axfc  +  Bu^  +  Dvfc 

yfc  »  Cxfc  +  k  >  0  (2.50) 
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(2.51) 


where  e  R*  and  w^  e  Rr  are  zero  mean  sequences  with 


Equation  (2.20)  applies,  E  is  nonsingular  and  D  e  Rnx*.  The  input 
sequence  u^  is  known,  and  the  initial  state  xq  is  uncorrelated  with 
and  w^. 

It  can  be  shown  that  the  optimal  observer  is  given  by 

Exk+1  -  Axr  +  Buk  +  K0^  -  Cxr)  .  (2.52) 

The  weighted  mean  square  reconstruction  error 

E{(*k  "  *k)Tw(*k  "  »  (2.53) 

where  W  >  0,  is  minimized  when  R  is  chosen  to  be 

K°  -  (AXCT  +  DS)(CXCT  +  R)-1  (2.54) 

and  X  is  found  by  solving  the  GARE: 

EXET  -  AXAT  +  DQDT  -  (AXCT  +  DS)(CXCT  +  R)-l(AXCT  +  DS)T  (2.55) 
or  equivalently 

EXE1  -  (A  -  DSR-lC)X(A  -  DSR_lC)T  +  DQDT  -  DSR_lSTDT 

-  (A  -  DSR-lC)XCT(CXCT+R)-lCX  (A  -  DSR_1C)T  (2.56) 

Equations  (2.55)  and  (2.56)  are  dual  to  (2.48)  and  (Z.47),  and  we  can 
solve  either  problem  using  the  same  technique.  We  shall  now  present  a 
technique  for  solving  the  GARE. 
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2.3.5  GARE  Solution  -  Continuous-Time  Case 


We  will  work  in  the  context  of  the  regulator  problem  for  conven¬ 
ience.  If  we  define 

M  E  0  A-BR_1ST  -BR_1BT 

y”[pJ;L“  oet  ;M“  sr_1st-ctqc  -(A-BR_1ST)T 


Then  we  can  associate  with  (2.29)  a  generalized  eigenproblem 
ALy  ■  My. 


(2.58) 


Theorem  2.3;  (Generalized  Real  Schur  Form)  let  L,  M  e  R2nx2n  and 
define  the  matrix  pencil  AL  -  M.  The  pencil  is  said  to  be  regular  if 
det(AL  -  M)  £  0.  There  exist  orthogonal  transformations  P  and  Z  such 


PT(AL  -  M)Z  -  APTLZ  -  PTMZ  «  AL  -  M 


(2.59) 


where  L  is  upper  triangular  and  ft  is  quasi-upper  triangular.  For  the 
lxldiagonal  blocks,  the  generalized  eigenvalues  are  real  A^  - 


(possibly  "infinite"  if  A a  -  0).  The  2x2  blocks  correspond  to  a 
finite  pair  of  complex  conjugate  eigenvalues.  Moreover,  these  eigen¬ 
values  can  be  arranged  in  any  desired  order. 

Proof :  [10] 

It  can  be  shown  for  L  and  M,  as  defined  in  (2.57),  that  if  A  is  a 
generalized  eigenvalue  of  AL  -  M,  then  -A  is  also  (Hamiltonian  property 
of  the  eigenvalues).  If  the  transformed  pencil  PT(AL  -  M)Z  is  parti¬ 


tioned  as 
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(2.60) 


where  L^,  e  R113™,  then  the  vectors  corresponding  to  the  first  n 
columns  of  Z  span  the  elgenspace  of  XL}}  -  [10].  Moreover,  we  can 

require  that  the  real  parts  of  the  generalized  eigenvalues  of  XL^  - 
Mu  be  negative.  If  we  partition  the  corresponding  Z  into  nxn  blocks 
as 


Z  - 


Z11 

Z21 


Z12 

Z22 


(2.61) 


we  can  state  the  following: 

Theorem  2.4:  With  respect  to  the  assumptions  made  for  this  generalized 
eigenvalue  formulation  with  L,  M  and  Z  as  defined  above, 

1)  Let  V  -  jj  Z  then  X  -  V21  Vu“l  (2.62) 

solves  (2.30)  with  X  »  X^  _>  0  ; 

2)  The  generalized  eigenvalues  of  XL^  -  are  the  closed- 
loop  eigenvalues  of  the  system  under  the  optimal  feedback 


u°(t)  -  -R"1 (S^  +  BTXE)x(t) 


Proof :  [12] 

The  above  method  depends  explicitly  on  R”1  and  can  encounter 
computational  difficulties  when  R  is  ill-conditioned  with  respect  to 
inversion.  If  this  is  the  case,  the  following  compression  technique 

due  to  Van  Dooren  [10]  can  be  employed.  Arrange  equations  (2.23) 

through  (2.25)  as 
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— 1 
w 

o 

o 

_ 1 

* 

A  0  B 

x 

0  ET  0 

- 

-ctqc  -at  -s 

p 

1 — 
o 

o 

o 

1 _ 

_  <S  _ 

_sT  bt  r_ 

_  u  _ 

(2.63) 


Determine  an  orthogonal  matrix  0  e  r(  2n+m)x( 2n4m)  such  that 


Ull 

1  - 

:  U12 

B 

-S 

m 

0 

0 

U21 

r°22 

R 

R 

(2.64) 


where  R  e  R™0™  and  Is  nonsingular.  Then  apply  U  to  the  pencil  corre¬ 
sponding  to  (2.63) 


(2.65) 


"e  0  0 

1 - 

« 

o 

< 

1 _ 

\ 

0  ET  0 

- 

-ctqc  -at  -s 

o 

o 

l _ 

sT  bt  r_ 

to  obtain  the  pencil 

e  0 


AUu 


0  et 


®ll 


A  0 
j_-CTQC  ~AT_ 


+  u12 


|^ST  B*j 


(2.66) 


It  can  be  shown  [10]  that  the  pencils  (2.66)  and  (2.58)  are 
equivalent.  Therefore,  the  pencil  (2.66)  can  be  used  to  solve  for  X 
without  having  to  Invert  R. 

2.3.6  GARE  Solution  -  Discrete -Time, Case 

As  in  the  previous  section,  we  work  in  the  context  of  the  regulator 
problem  for  convenience.  The  matrix  pencil  appropriate  for  this 
problem  can  be  formed  from  equations  (2.43)  through  (2.45)  as 


27 


Aj. 
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0  AT  0 
0  -Bt  0 


-CtQC  Et  -S 
ST  OR 


(2.67) 


Assuming  for  Che  moment  that  R  Is  Invertible,  we  can  derive  a  pencil 


equivalent  to  (2.67)  in  the  form  XL  -  M  as 


*  *1  P*  m 

E  BR-lBT  At-BR-1  St  0 

0  (A-BR“1St)t  -CtQC+SR_1 St  Et 


(2.68) 


where  L  and  M  are  defined  appropriately.  It  can  be  shown  that  the 
generalized  eigenvalues  of  (2.68)  have  the  symplectlc  property,  i.e., 
if  X  is  a  generalized  eigenvalue  of  (2.68),  then  is  also  [12]. 

We  can  now  solve  the  discrete-time  problem  in  a  manner  analogous 
to  the  continuous-time  case.  We  first  transform  the  pencil  (2.68)  to 
the  form  (2.60)  using  orthogonal  P  and  Z.  Note  if  R  Is  singular  or 
ill-conditioned  with  respect  to  inversion,  we  can  compress  the  pencil 
(2.67)  to  a  pencil  equivalent  to  (2.68)  using  the  technique  of 
Van  Dooren  [10].  Now,  we  will  require  that  Ixl  corresponding  to  the 


resulting  pencil  XL^  -  be  less  than  unity.  If  we  partition  Z  as 
in  (2.61),  we  have 

Theorem  2.5:  With  respect  to  the  assumptions  made  for  this  generalized 
eigenvalue  formulation  with  L,  M,  and  Z  as  defined  above. 


1)  Let  V 


E  0 
0  I 


X  -  V 


21  yll 


(2.69) 


solves  (2.48)  with  X  -  XT  J>  0; 


2)  The  generalized  eigenvalues  of  XL^  -  Mjj  are  the  closed- 
loop  eigenvalues  under  the  optimal  feedback 
u°k  -  -(R  +  BtXB)_1(St  +  BTXA)xk 

Proof:  f 12  J 

Note  that  in  this  discrete-time  problem,  neither  the  solution  of 
the  GARE  nor  the  optimal  feedback  depends  explicitly  on  R  ^  as  in  the 

continuous-time  problem.  The  special  case  R  *  0  is  called  "deadbeat 
control"  and  is  discussed  for  the  case  S  ■  0,  E  «  I,  utilizing  the 
above  methods  for  solution  in  [24]. 


CHAPTER  3 


ITERATIVE  REFINEMENT  BY  NEWTON'S  METHOD 

Numerical  implementation  of  the  GARE  solution  methods  of  Chapter  2 
Is  relatively  straightforward.  Proven  stable  algorithms  exist  that  are 
coded  into  reliable  FORTRAN  software  (see  Chapter  5).  However,  the 
GARE  may  be  ill-conditioned,  and  the  resulting  numerical  solution  may 
not  be  as  accurate  as  desired.  This  chapter  presents  an  iterative 
refinement  procedure  utilizing  Newton's  method.  This  procedure  is 
presented  for  the  continuous-  and  discrete-time  problems.  The  continu¬ 
ous-time  method  is  based  on  Klelnman  [13],  and  the  discrete-time  method 
is  based  on  Hewer  [14]. 

3.1  Iterative  Refinement  of  Continuous -Time  Solution 

The  following  result  is  due  to  Kleinman  [13]  for  the  controllable 
case,  which  was  extended  to  the  stablllzable  case  by  Sandell  [25]. 
This  result  is  for  the  ARE  and  is  presented  here  in  detail,  although 
from  a  slightly  different  approach,  to  provide  the  groundwork  for 
extension  to  the  GARE. 

Given  the  system  (2.18)  and  criterion  (2.19)  with  E*I,  OI  and 
S*0,  the  resulting  ARE  from  (2.27)  is 

ATX  +  XA  -  XBR_1BTX  +  CTC  -  Q.  (3.1) 

If  the  pair  [a,B]  is  stablllzable  and  the  pair  [C,A]  is  detectable, then 
there  exists  a  unique  solution  X-X^  >  0  to  (3.1),  such  that  the 


linear  feedback 
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u°  -  -R-1BTXx  -Kx  (3.2) 

stabilises  the  closed-loop  system  and  minimizes  the  criterion  (2.19). 
Furthermore,  it  can  be  shown  that 

o  T 

J(x0;u°)  -  min  J(x0;u)  -  x0Xx0.  (3.3) 

u 

To  derive  the  Newton  Iterative  procedure,  consider  at  the  k-th 
iteration  that  the  solution  X  of  (3.1)  is  of  the  form 

X  -  3^  +  6X  ,  (3.4) 

where  SX  is  "small",  i.e.,  can  neglect  second -order  terms  in  5X. 
Substituting  (3.4)  into  (3.1),  we  have 

0  -  ATXk+  XkA-XkBR-lBTXk-K:TC+6X(A-BR-1BTXlt)+(A-BR-lBTXk)T6X.  (3.5) 
Denote  the  solution  of  (3.5)  at  the  k-**4  iteration  by  6X^, and  let 

5Xk  -  Vi  -  *k  •  <3-6> 

Substituting  (3.6)  into  (3.5)  and  simplifying,  we  have 

0  -  (A-BV/Vl  +  +  °Tc  +  Cl“fcrt  •  <3-7) 

where 

,  (3.8) 

and  A-BK^+i  is  the  closed-loop-system  matrix. 

Equation  (3.7)  is  a  Lyapunov  equation  in  the  unknown  X^+i  and 


its  unique  non-negative  definite  solution  is  given  by 
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X]c+1  -  /“e(A"BKk+l>  C  (CTC  +  Kj+1RKk+1)e  (A-BKk+l)cdt ,  (3.9) 

which  Is  finite.  If  and  only  if,  the  closed  -  loop  -  system  matrix  is 
stable  (l.e.,  has  eigenvalues  with  negative  real  parts).  It  can  be 
shown  that 

X*  ■  *2  '  C  «“'8,C2,tt[<^^)T»(Kl-K2,-,Kl-K2,I(BTX1-8K2) 

-  (BTX1-RK2)T(K1-K2)]e(A_BK2)tdt  (3.10) 

or,  alternatively, 

T 

Xi  -  X2  -  /"  e(A"BKl)  t[(Kl-K2)TR(Ki-K2>  "  (K1-K2)T(BTX2-RK2) 

0 

-  (BTX2-RK2)T(K1-K2)]e(A“BKl)t  dt.  (3.11) 

We  can  now  state  and  prove  the  main  result  of  Kleinman  [13]  as 
extended  by  Sandell  [25]. 

Theorem  3.1:  Let  X^,  k»0,l,...,  be  the  unique  non-negative  definite 
solution  of  the  linear  algebraic  equation 

0  -  (A-Bl^)^  +  X^A-BI^)  +  CTC  +  (3.12) 

where,  recursively, 

1^  -  R-lBTXk_1  ,  k  -  1,2,...,  (3.13) 

and  Kq  is  chosen  such  that  the  matrix  (A-BKq)  is  stable.  Then, 

1)  0  <  X  <  Xfc+1  <  Xfe  <  ...  <  Xq 

2)  Xfc  m  X 

3)  in  the  vicinity  of  X,  IX^-XI  <  C2IXfc-XI2  , 
where  is  a  finite  constant. 
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Xq  -  Xj  -  /"  t(K0-K1)TR(K0-K1)e(A“BKl)t  dt  >  0  , 


so  Chat  Xj  <_  X(j .  In  addition,  we  have  by  (3.11) 


Xj  -  X  -  r  e(A“BKl)  t(K1-K)TR(Kl-K)e(A"BKl)t  dt  >  0. 


Hence,  Is  bounded  above  and  below  and,  therefore,  has  finite  norm. 


It  should  be  noted  that  another  iterative  scheme  is  possible  and 
has  been  used  in  practice  [26].  One  simply  solves  (3.5)  at  iteration  k 
for  and  then  takes 

Vi  ■  \  +  s\  •  <3-,6> 

It  has  been  shown  [27]  that 

iim^  53^-0  (3.17) 
and 

iim^  \  -  X  .  (3.18) 

This  procedure,  while  theoretically  equivalent  to  Theorem  3.1,  has  some 
computational  drawbacks.  Namely,  it  requires  more  machine  operations 
and  if  the  initial  guess  for  Xq  is  nonsymmetric,  then  X^+i  will  be 
non8ymmetric  at  each  iteration. 

Before  we  extend  Theorem  3.1  to  the  general  case,  we  first  want  to 
show  that  solving  the  GARE  (2.27),  which  results  from  the  generalized 
state  space  system  (2.18)  with  criterion  (2.19)  and  then  applying  the 
optimal  feedback  (2.31),  is  equivalent  to  solving  a  "standard"  regula¬ 
tor  problem  when  E  is  invertible. 

We  will  always  assume  that  E  is  invertible  since  when  E  is  singu¬ 
lar  the  GARE  (2.27)  does  not,  in  general,  have  a  solution.  However, 
even  though  E  is  assumed  nonsingular,  it  is  computationally  undesirable 
in  many  cases  to  invert  E  and  solve  the  equivalent  standard  problem.  A 
method  for  solving  the  generalized  regulator  problem  when  E  is  singular 
has  been  proposed  by  Cobb  [28]  and  involves  decomposing  the  problem 
into  a  standard  part  and  a  nonsingular  E  part. 


Theorem  3.2: 


Solving  the  generalized  regulator  problem  of  section 


P 


$ 

sfc 


I 


A 

vv 


u v, 


2.3.1  is  equivalent  to  solving  the  following  ''standard*’  problem 


x  “  Ax  +  Bu 
y  -  Cx, 

with  criterion 


(3.19) 


J  ■  i  /  (yTQy  +  2x^Su  +  u^Ru)dt 
*  n 


(3.20) 


where 


(3.21) 


(3.22) 


A  -  E-1  A 
B  -  E“lB. 

Proof:  The  ARE  corresponding  to  this  problem  is 

0  -  (A-BR~lST)TX  +  X(A-BR"lST)  -)5r-1BTX  +  ctqc  -  SR-lST 

-  (A-BR-lST)TE-TX  +  X  E_1 (A-BR_1ST)  -  XE-lBR-lBTE_TX 

+  CTQC  -  SR_1ST. 

The  optliaal  control  law  is 

u°  -  -R_l(BTX  +  ST)x 

-  -R-1 (BTE~TX  +  ST)x,  (3.23) 

which  results  in  the  closed -loop  system 

x  -  E“l(A-BR-lST  -BR-lBTE'TX)x.  (3.24) 

Now,  if  we  note  that 


X  -  etxe 


(3.25) 

them  (3.22)  is  equivalent  to  (2.27),  and  (3.24)  is  equivalent  to  (2.18) 
with  the  feedback  (2.31),  that  is., 
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Ex  -  (A-BR“lST  -BR"1BTXE)x. 


(3.26) 


Now  we  will  extend  Theorem  3.1  to  the  general  case. 

Theorem  3.3;  Let  X^,  k-0,1,...  be  the  unique  non-negative  definite 
solution  of  the  linear  algebraic  equation 

0  -  (A-BKk)TXlcE  +  ETXk(A-BK^)  +  CTQC  +  K^-  SK^  -  (SI^)1 


(3.27) 


where,  recursively. 


-  R"l(BTXk_1B  +  ST)  ,  k-1,2, . .. 


(3.28) 


and  where  Kq  is  chosen  such  that  the  matrix  E~*(A-BKq)  is  stable*  Then 

n  o<x<xk+1<xk<...<x0 

2 )  %  imj^  ■  A 

3)  In  the  vicinity  of  X,  lXfc+1  -  XI  <  ^IX,^  -  XI2 

where  X  solves  the  GARE  (2.27),  and  Cj  19  a  finite  constant. 

Proof;  1)  We  employ  the  results  of  Theorem  3.2  to  convert  the  problem; 


from  (3.25)  let 


xk  -  e-\e-i 


(3.29) 


and  substitute  In  (3.27)  and  (3.28)  to  obtain 


0  -  [ E”1 ( A-BK^ ) ] TX^  +  X^E-1 (A-BK^)  +  CTQC  +  K^RK^ 


-  -  <SKfc)T, 


(3.30) 


where 


\  -  R-1{(E-lB)TXk._1  +  ST} 


(3.31) 
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Xj.  -  /%tE"l<A"B'k)lTt(CTqC  +  kJr^-  -  (S^)1) 

eE-‘(*-B\)V 

when  E-1(A-BK.  )  Is  stable.  It  can  be  shown  that 


(3.32) 


Xi  -  X2  -  r  JE’1(A-B*2>1  t[(K1-K2)TR(K1-K2) 

0 

-  (K1-K2)T{(E-1B)TX1  -  rk2  +  sT} 

-  { (E“1B)TXj  -  RK2  +  ST}T(Kl-K2)]eE  l(A_B^2>tdt 


(3.33) 


or,  alternatively. 


-  Xj ,  -  /*  e  f  E~ 1 ( A_BEl } J  ^ ( ( Ki -K2 ) TR( Ki -K2 ) 

0 

-  (K1-K2)T{(E-1B)TX2  -  rk2  +  sT{ 

-  {(E-lB)TX2  -  RK2  +  ST}T(Kl-K2)]eE  1(A‘®Kl),:dt  . 

(3.34) 

Now  let  Xg  satisfy  (3.30)  for  a  Kg  chosen  such  that  E-1(A-BK0)  Is 
stable.  We  remark,  here  that  by  Theorem  3.2,  this  Is  equivalent  to 

stabilizing  E-^A-BKq).  Set  Kx  -  R_1  ((E_1B)TX0  +  ST),and  let  XL  be  Che 
associated  solution  to  (3.32).  Using  (3.33),  we  obtain 


-  *i  -  r  • 


[E-1  (A-BK, )  ]Tt  [  ( Jo.Si  ,IR( Ro_Rl)  )eE-‘  (A-BK,  )t 


(3.35) 


so  that  Xj<  Xo.  In  addition,  we  have  by  (3.34) 
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Xi  -  X  -  /“  e[E‘l(A-BKl>]  t[(K1-K)TR(Ki-K)]eE"1(A_B^l)tdt  >  0. 

(3.36) 

Hence.  Xj  Is  bounded  above  and  below  and,  therefore,  has  finite  norm. 

Thus,  E“*(A-BKi)  is  stable  so  Xj  satisfies  (3.30)  with  k«l.  Repeating 
the  above  argument  for  k-2,3,...  yields 

0  <  *  <  \+1  <  \  <  •••  <  *o 


Since  by  (3.29)  E  is  a  congruence  transformation  on  X  to  yield  X,  the 
desired  result  is  implied. 

2)  Taking  the  limit  of  (3.27)  as  k-*»,  we  obtain  after  some 
manipulation 

0  -  (A-BR“lST)TX  E  +  ETX  (A-BR-^S1)  -  etx  br-1btx  E 

OD  OB  OB  OB 

+  CTQC  -  SR^S1  .  (3.37) 

Since  X  is  the  unique  non-negative  definite  solution  to  (3.37),  X  »  X. 

OB 

3)  Set  iq  -  R-1  [(E^B)1^  +  ST]  and  X2  -  R-l[(E_lB)TX  +  ST] 
in  (3.34)  to  find 

Xk+1  -  X  -  /"e^E  C  (Xk-X)E-1BR"1BTET 

(Xk-X)eE~1(A_B^k+l)t  dt.  (3.38) 

Substitute  (3.25)  into  (3.38)  to  obtain 
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ET(Xk_fl-X)E  -  /"e^E  1  (A_BKk+l *  1  tET(Xk-X)BR-lBT 

(Xjt~X)EeE”1  ^A_BISc+l^t  dt.  (3.39) 

-T  i 

Pre-multiply  (3.39)  by  E  ,  post-multiply  by  E  ,  and  take  the  norm  to 
find 

IX^  -  XI  <  /"leE  1(A_BKk+l)tl2dtlE-1l2lEI2IBR-1BTIIXk-XI2. 

(3.40) 

It  can  be  shown  that,  uniformly  In  k, 

C  leE  BKk+l^Cl2dt  <  constant  ■  . 

o  — 

Let  C2  -  C1k2(E)IBR“1BtI,  and  the  proof  is  complete. 

tfe  remark  that  the  proof  to  part  3  Indicates  that  the  square  of 
the  condition  of  E  with  respect  to  inversion,  ic2(E),  Influences  the 
convergence  rate  of  the  algorithm  oultlplicatlvely. 

3.2  Iterative  Refinement  of  Discrete-Time  Solution 

An  iterative  scheme  is  possible  for  the  discrete-time  formulation, 
and  it  was  first  reported  by  Hewer  [14J.  The  result  is  for  the  dis¬ 
crete-time  ARE,  and  the  proof  given  here  is  similar  to  the  continuous- 
time  case. 

Given  the  system  (2.40)  with  criterion  (2.41)  with  E-I,,Q-I,  and 
S»0,  the  resulting  ARE  from  (2.48)  is 

X  -  ATXA  -  ATXB(R+BTXB)_1BTXA  +  CTC  .  (3.41) 

If  the  pair  [a,b]  is  stabilizable  and  the  pair  [C,a]  is  detectable, then 
there  exists  a  unique  solution  X  ■  X^  >  0  to  (3.41)  such  that  the 
linear  feedback 

40 


V  '  f  o'.  >T.  •*. 
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u°  -  -(R+BTXB)-lBTXAxk  -Kxk  (3.42) 

stabilises  the  closed-loop  system  and  minimizes  the  criterion  (2.41). 

One  can  follow  a  first-order  perturbation-type  derivation  to 
arrive  at  the  Lyapunov  equation  that  forms  the  basis  of  the  following 
theorem  in  a  manner  similar  to  that  of  the  previous  section  for  the 
continuous-time  problem.  The  following  result  is  due  to  Hewer  [14]. 
Theorem  3.4:  Let  X^,  k-0,1,...,  be  the  unique  non-negative  definite 
solution  of  the  linear  algebraic  equation 

\  -  (A-BKk)TXk(A-BKk)  +  +  CTC,  (3.43) 

where,  recursively, 

\  -  (R+B^jBr^X^A,  k-1,2 .  (3.44) 

and  Kg  is  chosen  such  that  the  closed -loop -system  matrix  (A-BR0)  is 
stable  (i.e.,  has  eigenvalues  whose  magnitudes  are  less  than  unity). 
Then 

D  °<x<  Xk+i  <  \<  ...  1*0 

2>  u»k-  \  ’  x 

3)  in  the  vicinity  of  X,  IXk+1  -  XI  <  C2IXk  -  XI2. 
where  C2  is  a  finite  constant. 

Proof :  1)  If  (A-BKfc)  is  stable,  the  unique  non-negative  definite 

solution  Xfc  of  (3.43)  may  be  written  as 
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(3.45) 


\  -  „I0  [(A-BKk)T]N(KjRKk+CTC)(A-BKk)N  . 

It  can  be  shown  chat  for  (A-BK^ )  and  (A-BK2)  stable  then 

Xx  -  X2  -  (A-BK1)T(X1-X2)(A-BK1)  +  (K1-K2)T(R-«TX2B)(K1-K2) 

+  [(R+BTX2B)K2-BTX2A]T(K1-K2)  +  (K1-K2)T[(R+BTX2B)K2 
-BTX2A]  (3.46) 

-  jjIq  [(A-BK1)T]N((K1-K2)T(R+BTX2B)(Kl-K2)+[(R+BTX2B)K2 

-BTX2A]-KX1-K2)T[  (R+BTX2B)K2-BTX2A])  (A-BKj  )N,  (3.47) 

or,  alternatively, 

Xi  -  X2  -  (A-BK2)T(X1-X2)(A-BK2)  +  (K1-K2)T(R-WTX1B)(K1-K2) 

+  [  (R+BTXjB)K2-BTX1A]T(K1  -K2  )+(  XL  -K2  )T[  (R+B^!  B)K2 
-BTXiA]  (3.48) 

-  nE0  [(A-BK2)T]N((Ki-K2)T(R+BTX1B)(K1-K2)+[(R+BTX1B)K2 

-BTX1A]T(K1-K2)  +  (K1-K2)T[(R+BTXlB)K2-BTX1A]) 

(A-BK2)N  .  (3.49) 

T  1 

Let  Xg  satisfy  (3.43)  for  the  chosen  Kg.  Now  set  K*  ■  (R+B  X0B)~A 

BTXgA  and  let  X|  be  the  associated  solution  to  (3.43).  Using  (3.49)  we 
obtain 

Xg  -  Xi  -  N|0  [(A-BK1)T]N[(K0-K1)T(R+BTX0B)(K0i-K1)](A-BK1)N>0, 
so  that  Xi<Xg.  In  addition,  we  have  by  (3.47) 

XL  -  X  -  JL  [(A-BKl)T]N[(Kl-K)T(R+BTXB)(K1-K)](A-BKl)N  >  0. 


Hence,  Xj  Is  bounded  above  and  below  and,  therefore,  has  finite  norm. 
Thus,  (A-BKj)  Is  stable  so  X}  satisfies  (3.43)  with  k»l.  Repeating  the 
above  argument  for  k>2,3,...  yields  the  desired  result. 

2)  Talcing  the  limit  of  (3.43)  as  k-n»,  we  obtain 

X  *  ATX  A  ■  ATX  B(R+BTX  B]flBTX  A  +  CTC.  (3.50) 

00  00  00  00  00 

Since  X  is  the  unique  non-negative  definite  solution  of  (3.50),  Xa,**X. 

3)  Set  Ki  -  (R-WTXkB)~lBTXkA  and  K2  -  (R+BTXB)~lBTXA  in  (3.47), 
and  take  the  norm  to  find 

IXfc+1-  XI  <  IXk-  XI2lB(R+BTXB)~1BTlN|0  l(A-BIk+l)*fll2. 

Since  (A-BK^+i)  is  stable,  it  can  be  shown  that,  uniformly  in  k, 

N|0  I  (A-BK^j)1*1!2  <  constant  -  Cx 

Let  C2  ■  C1IB(R+BTXB)“lBT|> and  the  proof  is  complete. 

It  should  be  noted  that  another  iterative  scheme  is  possible  here, 
as  in  the  continuous  case,  where  one  solves  for  a  correction  term  to 
the  solution  at  each  iteration.  This  procedure  faces  the  same  computa¬ 
tional  drawbacks  as  in  the  continuous  case  and,  therefore,  will  not  be 
discussed  further. 

Before  we  extend  Theorem  3.4  to  the  general  case,  we  first  want  to 
show  that  solving  the  GARE  (2.48),  which  results  from  the  generalized 
state  space  system  (2.40)  with  criterion  (2.41)  and  then  applying  the 


optimal  feedback  (2.49),  is  equivalent  to  solving  a  "standard"  regula¬ 
tor  problem  when  E  is  invertible. 


Theorem  3.5:  Solving  the  generalized  regulator  problem  of  section 
2.3.3  la  equivalent  to  solving  the  following  ‘’standard"  problem 


*k+l  "  **k  +  K 

yk  *  CXfc  (3.51) 

with  criterion 

J  "  7  klo  (yk<*k+  2xksV  uJRuk)  (3-52) 

where 

A  -  E“lA 

B  -  E~lB  .  (3.53) 

Proof :  The  ARE  corresponding  to  this  problem  is 

X  -  ATXA  -  (ATra+S)(R+BTS)“1(ATS+S)T  +  ctqc 

-  ate-txe_1a  -  ( ae_txe_ 1 B+S ) ( r+bte_txe~ 1 B ) - 1 ( ate_txe~ 1 B+S ) t 

+  CT0C  .  (3.54) 

The  optimal  control  law  is 

u°  -  -(R+BTXB)-l(ATXB+S)Txk 

-  -(R+BTE"TXE-lB)-l(ATE’TxE-lB+S)Txk  (3.55) 

which  results  in  the  closed-loop  system 

xk+1  -  E"l[A-B(R+BTE"TXE'lB)”l(ATE"TXE-lB+S)T]xk  .  (3.56) 


Now,  if  we  note  that 


irWlVi ■ '.’VW ^  ^  ».  V  *- 
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then  (3.54)  Is  equivalent  to  (2.48),  and  (3.56)  is  equivalent  to  (2.40) 
with  the  feedback  (2.49);  that  is, 


Exk+1  -  [A-B(R+BTXB)-1(ATXB+S)I]xk. 


Now  we  will  extend  Theorem  3.4  to  the  general  case. 


(3.58) 


Theorem  3.6:  Let  X^,  k-0,1,...  be  the  unique  non-negative  definite 
solution  of  the  linear  algebraic  equation 

E^E  -  (A-BKk)TX|c(A-BKIc)  +  +  CTQ C 


-  -  (SK^) 


where ,  recursively , 


-  (R+BTXklB)-1(BTXk_1A+ST),  k-1,2,.... 


(3.59) 


(3.60) 


and  Kg  is  chosen  such  that  the  closed  loop  system  matrix  E-I(A-BK0)  is 


stable.  Then 


1)  o  <  x  <  xk+1  <  xk  £  ...  <  * 

2)  -  X 

3)  in  the  vicinity  of  X,  |Xk+1  -  XI  <_  CalX^  XI2 

where  X  solves  the  GARE  (2.48),  and  C2  is  a  finite  constant. 

Proof :  1)  We  employ  the  results  of  Theorem  3.5  to  convert  the 


problem;  from  (3.57)  let 

Xk  -  e-\e-i 

and  substitute  in  (3.59)  and  (3.60)  to  obtain 


(3.61) 


W3CWW’ 
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\  “  [E~*  (A-BKyt)]TX^E-^  (A-BK^)  +  +  CTQC 


where 


-  SK^  -  (SK^)1 


Ky.  -  (R+BTE~TXk_1E-1B)-1(BTE'TXlc_1E-1A+ST) 


0.6: 


\  "  1S0  [(E-I(A-BKk))T]N[KjRKk4CTQC-SKk-(SKk)T][E-1(A-BKk)]* 


(3.6^ 


when  E-1  (A-BK^,)  is  stable.  It  can  be  shown  that 


Xi  -  *2  -  n?0[(e-1(a-bk1))t]n[(k1-k2)t(r+bte'tx2e-1b)(k1-k2) 


*  v.r 


+  (  (R+BTE*'TX2E-lB)K2-BTE~TX2E-1A-ST)T(K1-K2) 
+  (R-rt^^^E-Uj^-B^'^E-^-S1)  ] 

[E-1 (A-BK1)]N 


(3.65 


or,  alternatively 


X!  -  Xj  -  J0  [(E-l(A-BK2))T]N[(K1-K2)T(R+BTE“TX1E-lB)(Kl-K2) 


+  ((R+BTE*TXlE-1B)K2  -  B^^Xj E-1  A  -  ST)T(K1-K2) 

I 

+  (Kx  -itj  )T(  (R+BTE~TXl  E-1  BJ^-B^'^E-^-S1)  ] 


[E-1 (A-BK2)]N. 


(3.66 


NWC  TP  6521 


Now  let  Xg  satisfy  (3.62)  for  a  Kg  chosen  such  that  E-i(A  -  BKq)  is 
stable.  Let  X}  be  the  associated  solution  to  (3.64).  Using  (3.66)  we 


obtain 


Xq  -  Xj  -  N|0  [(e-1(a-bk1))t]n[(Ko-Ki)T(R+bVtx0e-1b)(k1-k2)] 

[E“l(A-BK!)]N  >  0 

so  that  Xj  £  Xq •  In  addition,  we  have  by  (3.65) 

Xi  -  X  -  N|0  [(E-1  (A-BKi  ))T]N[  (X1-K)T(R+bVTXE-1B)(K1-K)  ] 

[E-^A-BKi)]11  >  0. 

Hence,  X1  is  bounded  above  and  below  and,  therefore,  has  finite  norm. 

Thus,  E-* (A-BK^ )  is  stable  so  X^  satisfies  (3.62)  with  k-l.  Repeating 
the  above  argument  for  k~2,3,...  yields 

0  £  X  £  Xk+1  £1^  £  ....  <  X0  . 

Since  by  (3.61)  E  is  a  congruence  transformation  on  X  to  yield  X,  the 
desired  result  is  implied. 

2)  Taking  the  limit  of  (3.59)  as  k+»  we  obtain 

ETX  E  ■  ATX  A  -(ATX  B+S)(R+BTX  B)_1(ATX  B+S)T  +  CTQC  .  (3.67) 

00  00  00  00  A 

I 

Since  X  is  the  unique  non-negative  definite  solution  to  (3.67),  X^^X. 

3)  It  is  easily  shown  from  (3.65)  that 


'O.v'v'Oa'  OO 


ET(XrX2)E  -  N|0  [(A-BK1)T]N[(K1-K2)T(R+BTX2B)(K1-K2) 

+  ((R+BTX2B)K2  -btx2a-st)t(k1-k2) 

+  (K1-K2)T((R+BTX2B)K2  -BTX2A-ST)][A-BK1]N. 

(3 

Set  -  (R+BTXkB)-1(BTXkA+ST)*and  K2  -  (R+BTXB)_1 (BTXA+ST) 

in  (3.68),  and  take  norms  to  find 

IXfc+1  -  XI  <  IXk-XI2  lE-1l2IB(R-fBTXB)-lBTlN|0l(A-BKk+1)N+1 
Since  (A-BKfc+i)  is  stable,  it  can  be  shown  that,  uniformly  in  k, 

^  l(A-BKk+1)l*fll2  <  constant  -  Cj . 

Let  C2  «  Cj  IE"1I2IB(R+B1'XB)“1BT|,  and  the  proof  is  complete. 
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CHAPTER  4 

CONDITIONING  OF  THE  ALGEBRAIC  RICCATI  PROBLEM 

A 8  stated  in  Chapter  2, It  Is  desirable  to  associate  a  measure  with 
a  computing  problem  which  reflects  the  overall  sensitivity  of  the 
solution  to  changes  in  the  data.  It  is  the  intent  of  this  chapter  to 
derive  such  a  measure  (i.e.,  condition  number)  for  the  problem  of 
computing  a  solution  to  the  GARE.  To  be  most  useful,  the  condition 
number  must  be  easily  computable.  It  should  not  require  additional 
computations  on  the  same  order  of  computing  the  solution  itself. 

We  will  first  state  the  desired  form  for  the  condition  number  and 
review  previous  work  on  condition  estimates  for  the  Rlccatl  problem. 
Then  first  order  perturbation  analysis  will  be  employed  to  derive  new 
condition  estimates  for  the  continuous-  and  discrete-time  Rlccatl 
equations.  Ways  of  employing  balancing  to  improve  the  numerical 
accuracy  of  the  calculated  solution  will  be  discussed.  An  efficient 
method  for  incorporating  a  change  of  model  coordinates  is  presented. 

4.1  Previous  Work 

We  seek  a  relative  condition  estimate  for  the  Rlccatl  problem  in 
the  following  form: 

,  (4.1) 

where  X  is  the  exact  solution  to  the  GARE,  X  is  the  computed  solution, 
C  is  a  constant  possibly  depending  on  the  size  of  the  problem,  k(X)  is 
the  desired  condition  number,  and  em  is  the  machine  epsilon  (preci¬ 
sion).  Machine  epsilon  is  defined  as  the  smallest  positive 
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number  that  can  be  added  to  1.0  such  that  the  machine  recognizes  that 
the  sum  is  different  from  1.0.  In  this  form,  if  k(X)  is  0(10N), 

then  the  computed  solution  can  be  expected  to  differ  from  the  true 
solution  by  N  significant  digits. 

He  should  note  here  that  this  type  of  indicator  of  solution 
accuracy  is  desired  because  the  more  traditional  method  of  examining 
the  size  of  the  residual  is  not  always  reliable.  The  residual  is  the 
remainder  quantity  obtained  when  the  computed  solution  is  substituted 
in  the  original  problem.  That  is,  for  the  GARE  (2.30)  the  residual  is 

residual  (A-BR-1 ST)TXE  +  ETX(A-BR_lST)  -  ETXBR“lBTXE 

+  CTQC  -  SR”1 .  (4.2) 

Stewart  [15]  explores  in  detail  the  behavior  of  the.  residual  in  the 
linear  equations  problem  Ax-b  and  cites  examples  where  the  residual 

b  -  Ax  is  not  a  reliable  indicator  of  solution  (x)  accuracy.  There  is 
no  reason  to  expect  that  the  more  complex  problem  of  the  GARE  produces 
residuals  that  are  better  behaved,  even  though  counter-examples  are 
more  difficult  to  construct. 

We  look  for  k(X)  of  the  form: 

<(X)  -  f (A,B,C,E,Q,R,S)  .  (4.3) 

That  is,  k(X)  should  be  a  function  of  the  matrices  Involved  in  the 
GARE.  We  want  the  norms  involved  in  (4.1)  or  in  the  computation  of 
ic(X)  to  be  the  1,  •  or  Frobenius  norm  versus  the  2  norm  to  minimize  the 
number  of  calculations  involved. 


tassaftaa  xxxjlv&l  *&&&&>  &&j&ml 
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There  have  been  several  attempts  to  derive  bounds  for  the  solution 
to  the  ARE  [  29 ] — [  34 ] ,  none  of  which  are  of  the  form  (4.1).  The  follow¬ 
ing  result  by  Byers  [20]  is  for  the  continuous-time  ARE  (2.27): 


•cB(X) 


lCT0Cl  +  2lA||XI-HBR~1BTHXa2 
IX«SEP[A^,-AC] 


(4.4) 


where 


A  -  A  -  BR"1BTX 
c 


(4.5) 


SEP(F.G) 


a  pi-1 


JPF  -  GPI  . 


(4.6) 


Byers  arrives  at  this  bound  via  a  first -order -perturbation  approxima¬ 
tion  to  (2.27).  A  notable  feature  of  this  condition  estimate  is  that 
it  includes  the  effect  of  the  separation  of  the  closed -loop  spectrum 
(4.6)  on  the  condition  of  the  problem.  For  a  detailed  definition  of 
SEP(F,G)  and  a  discussion  of  its  properties  see  Stewart  [35], [36]. 
Some  speculation  [6],  [ 20 ]  and  empirical  results  suggest  that  a  small 
separation  of  the  closed -loop  spectrum  causes  loss  of  numerical  accur¬ 
acy  in  the  computed  solution  (see  Chapter  5).  A  potential  drawback  to 
Byers'  result  (4.4)  is  the  appearance  of  I XI 2  in  the  numerator  and  I XI 
in  the  denominator.  This  would  suggest  ill -conditioning  of  the  problem 
when  I XI  is  large  or  small,  which  is  not  necessarily  true. 

Martensson  [5]  showed  that  V^,  which  must  be  inverted  to  form  the 
Rlccati  (2.27)  solution,  i^  singular  if  the  model  is  unstabilizable  for 
the  case  when  eigenvectors  are  used  as  the  basis  for  the  stable  eigen- 
space  in  the  solution  process.  Laub  [7],  Pappas,  et.al.  [8],  and 
Emami-Naeinl  [9]  prove  analogous  results  when  Schur  vectors  are  used  as 
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the  basis  for  the  stable  eigenspace .  This  fact  and  the  fact  that  V}  j 
oust  be  inverted  (that  is,  a  linear  system  of  the  form  XVj^  «  V2i  must 
be  solved  for  X)  has  prompted  speculation  [6],  [7]  that  the  condition 
of  Vu  with  respect  to  inversion,  k(Vj^),  might  be  a  good  indicator  of 
the  condition  of  the  Riccati  problem  (2.27).  Indeed,  empirical  results 
show  that  when  Vj^  is  ill-conditioned  with  respect  to  inversion,  say 
k(Vi})  is  0(10*);  then,  in  general,  N  digits  of  accuracy  are  lost  in 
the  solution  X.  Since  for  the  GARE  solution,  inversion  of  involves 
the  Inversion  of  E;  then  <(V^)  can  be  expected  to  reflect  any  ill- 
effects  that  near  singularity  of  E  might  have  on  the  solution.  How¬ 
ever,  examples  do  exist  (see  Chapter  5)  where  the  Riccati  problem  is 
ill-conditioned,  but  this  fact  is  not  indicated  by  the  magnitude  of 
k(Vu). 

Another  indicator  of  numerical  accuracy  in  the  solution  of  the  ARE 
has  been  suggested  by  Paige  and  Van  Loan  [37].  Their  reasoning  is 
based  on  the  following  result: 

Theorem  4.1:  Let  the  unitary  matrix  E  c  Cnxn  be  partitioned  in  the 
form 

w.|wn  wi 

1*2 1  *2 

where  V\\  e  Cmxm  with  m  £  Then  there  are  unitary  matrices  U  - 

dlag  (UlfU2)  and  V  -  diag(V1,V2);  Vl  .Vj,  e  C®» »;  U2,V*  e 
such  that 

n  r  -A  0 

lAfV  -  A  r  0  (4.7) 

0  0  1 

52 
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where 

T  -  dlag  (Yi ,Y2 . • • • .Ym)  >  0 

A  m  dlag  (6i  *60  i  •  •  •  »6  )  )  0 

o  — 

Yi  <  Y2  <  •  •  •  <  Y„  <  1  •  (4.8) 

Proof:  [38] 

The  use  of  this  theorem  in  [37]  suggests  that  if  we  perform  this 
special  singular  value  decomposition  on  the  orthogonal  Z  found  In  our 
reduction  of  the  generalized  elgenproblem  (2.59)  we  have 
Theorem  4.2:  Given  Zjj ,  Z2j  as  found  in  (2.59),  (2.61) 

1)  Then  there  exists  orthogonal  U  and  V  e  R11*11  such  that 

UTZUV  -  r  -  dlag  (Y1.Y2 . Y„)  5 

T 

U  Z21V  ■  A  ■  dlag  (6j ,62 » ••  •  »6n)  » 

Yt2  +  i*  -  1;  0<Tl  <Y2<---<Yn<l  •  (4.9) 

6  _ 

2)  If  T  Is  nonsingular,  X  ■  U  diag(— , ...»  — )  U  solves  the 

1 1  Yn 

ARE  (2.27). 

Proof:  [37] 

The  size  of  Yi  Is  suggested  in  [37]  as  being  an  indicator  of 

potential  numerical  difficulties  in  the  solution  of  the  ARE  because 

.»  1/2 

IXI2  -  ■  U-Yi2)  /Yl  *  cotan  Yl  (4.10) 

and  if  X  Is  computed  on  a  computer  having  a  machine  precision  e,  then 
rounding  erroru  of  order  zh\  can  be  expected  to  contaminate  the 
result.  Therefore,  1/yi  may  indicate  the  condition  of  the  ARE. 
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We  can  relate  1/yi  to  ^(Zn)  (k(Zji)  -  K(Vn)  for  the  ARE  (2.27)), 
which  was  suggested  previously  as  a  condition  measure.  Recall, 

««n>  -  ,zn'2  •zii-l*2 

•  ”*<°<2n>>  sTOFTZITTT  (4-u> 

where  o(2j^)  denotes  the  set  of  singular  values  of  2^  [15].  There¬ 
fore, 

<(ZU)  -  ^  <  L.  .  (4.12) 

So  we  see  the  two  estimates  are  of  the  same  order  unless  all  of  the 
singular  values  of  Zji  are  small. 

4.2  Flrat-Order-Perturbatlon  Analysis  of  the  GARE 

To  facilitate  our  analysis  of  the  continuous- time  Rlccatl  problem 
in  this  section,  we  start  with  the  GARE  (2.30)  and  make  the  following 
definitions: 

D  A-BR-1 ST 

G  BR~*  BT 

H  CTQC-SR”lST  .  (4.12) 

The  GARE  can  now  be  written  as 

0  -  DTXE  +  ETXD  -  ETXGXE  +  H  F  (X)  .  (4.13) 

c 

Suppose  that  there  is  a  small  perturbation  in  F  (X),  say  5F  ;  we  would 

c  c 

like  to  know  the  effect  of  this  perturbation  on  the  solution  X. 

Theorem  4,3:  For  consistent  matrix  norms,  and  given  a  small  perturba¬ 
tion  5F  to  (4.13),  we  have  for  the  resultant  change  in  the  solution  6X 
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K 


; 


« 

■ 


E 


V. 

»v 


£ 


& 


K 


16X1  >  __ 


I6F  I 

c 


IE  11A  1+lEllA  I 

c  c 


(4.14) 


and 


16X1  < 


k(E)ic(E  )I6F  I 
c 


I E 1 1 ET I SEP[ Ac  E_1 , -( Ac  E- 1  )T ] 


(4.15) 


where 


k(E)  :■  |EI1E”*|  ■  condition  of  E  with  respect  to  inversion 


D  -  GXE  ■  closed -loop -system  matrix. 


Proof :  Starting  with  (4.13),  define 


F  (X-6X)  F  (X)  +  6F 
c  c  c 


(4.16) 


DT(X-6X)E  +  ET(X-6X)D  -  ET(X-6X)G(X-6X)E  +  H. 


Neglecting  second -order  terms  in  6X, 


Fc(X-6X)  -  DTXE  +  ETXD  -  ETXGXE  +  H  -  (D-GXE)T5XE 


-  E  6X(D-GXE) 


-  Fc(X)  -  (D-GXE)T6XE  -  ET5X(D-GXE) 


Therefore, 


-6F  -  A  T5XE  +  ET6XA 

c  c  c 


(4.17) 


Taking  norms  in  (4. 17),  we  have 


I6F  I  <  HAT  U5XIIE1  +  9ET|I6X9IA  I 
c  —  c  c 


-T 


which  yields  (4.14).  We  pre-multiply  (4.17)  by  E  and  post-multiply 


7-1 


by  E  to  obtain 
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-E~T5F  E-1  -  (A  E_1)T6X  +  5XA  E"1  . 
c  c  c 

Taking  norms  we  have 

I ZA  E"1 +(A  E— 1 )TZI 

IE  eF  E-1!  -  - £ - £ _  16X1 

C  IZf 

6  X 

where  Z  • -  .  Therefore, 

16X1 


lE'TUE-1il6Fcl  >  JjJ.j 


I  PA  E_1+(A 
c  c 


e-1)tpi 


16X1 


I  PI 


-  I6XISEP[AcE“1,  -(ace-1)t] 

which  yields  (4.15). 

Theorem  4,4;  Given  Che  GARE  (2.30),  a  relative  condition  number  for  the 
solution  X  is 


kA  (X) 
c 


k (E)k(ET ) I CTQC-SR- 1 ST| 

ixiieiieti  sep[ace_1  .-(i  r1)1] 


(4.18) 


where  A  .  (A-8I"1ST-Bi“1BTXB). 
c 

Proof;  To  facilitate  the  analysis,  we  consider  perturbations  in  FC(X) 

on  the  order  of  a  perturbation  in  the  constant  term,  i.e., 

I6F  I  <  IHI  •  e  ,  (4.19) 

c  —  m 

so  that  from  Theorem  4.3  we  obtain 


16X1 

T XT’ 


<(E)<(E  )|H| 


IXIIEIIETI  SEP[AcE-1 ,-(AcE-1'T 


e  . 

m 


(4.20) 


From  (4.20)  we  can  make  the  following  definition: 


:*  relative  condition  number  for  the  continuous-time 
GARE. 

We  remark  here  that  the  condition  number  of  Theorem  4.4  has  the  advan¬ 
tage  that  it  includes  the  effects  of  ill -conditioning  of  E  with  respect 
to  inversion,  the  separation  of  the  closed-loop  spectrum  anid  the 
ill -conditioning  of  R  with  respect  to  inversion  when  S  is  non-zero.  It 
does  have  the  disadvantage  of  having  I XI  in  the  denominator,  as  did 
Byers'  result.  Clearly,  the  form  of  the  bound  on  the  perturbation, 
6FC,  dictates  the  form  of  the  resulting  condition  number.  One  could 
consider  perturbations  to  all  matrices  Involved  in  the  problem  with  a 
corresponding  increase  in  the  complexity  of  the  resulting  condition 
number.  However,  at  this  point  it  is  not  clear  that  this  is  necessary, 
and  it  is  counter  productive  to  our  goal  of  a  simple,  easily  computable 
measure.  The  proper  form  of  the  bound  on  (Fc  is  a  topic  of  continu¬ 
ing  research. 

We  now  turn  our  attention  to  the  discrete-time  Riccatl  problem. 
To  facilitate  our  analysis,  we  start  with  the  GARE  (2.48)  and  make  the 
following  definition: 

F.(X)  ATXA  -  ETXE  -  (ATmS)(R4BTXB)-1  (ATXB+S)T 
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As  la  the  continuous-tine  esse,  suppose  there  is  a  snail  perturbation 


in  F^(X),  say  6F4,  we  would  like  to  know  the  effect  of  this  perturba¬ 


tion  on  the  solution  X. 


Theorem  4.5;  For  consistent  matrix  norms,  and  given  a  small  perturba¬ 


tion  6F.  to  (4.22),  we  have  for  the  .resultant  change  in  the  solution  6X 
d 


16X1  > 


IE4IIEI  +  IA.  1 I A. I 
a  a 


(4.23) 


16X1  < 


»6Fdl 


-  ,AdT,«Ad. 
<(E)k(ET)  d  d 


(4.24) 


IEIIET|  >  ^EME^IA^nAjI 


(4.25) 


where 


k(E)  :■  IEIIE-1I  «  condition  of  E  with  respect  tQ  inversion 


(A-BK)  -  closed -loop -system  matrix 


K  (R+BTXB)"1 (ATXB+S)T  . 


Proof:  Starting  with  (4.22)  define 


Fd(X+6X)  Fd(X)  +  6Fd 


-  AT(X+6X)A  -  ET(X+6X)E  +  CTQC 


-  [ A1 (X+6X)B+S] [ R+B 1 (X+6X)B ] -l [ AT(X+6X)B+S ] 


After  expanding,  neglecting  second-  or  higher-order  terms  in  5X,  and 


some  algebraic  manipulation  we  obtain 


$ 
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Fd(X+6X)  -  Fd(X)  -  ET6XE  +  (A-BK)TSX(A-BK) 


Therefore, 


6F.  -  -ET6XE  +  AjTSXA.  . 
a  a  a 


Taking  norms  in  (4.26)  we  have 


(4.26) 


l6Fdl  <  |ET| 16X1 IEI  +  IAdTfl6XllAdl 

which  yields  (4.23).  We  pre-multiply  (4.26)  by  E  and  post-uultlply 
by  E_1  to  obtain 

E""T6F,E“*  -  -«X  +  (AjE“1)T6XAjE-1. 
a  ad 


6X  -  (A.E”1  )T6XA  E-1  -  E^FjE-1  . 
a  a  a 

Taking  norms  yields  (4.24). 

Theorem  4.6:  Given  the  GARE  ( 2 . 48),  a  relative  condition  number  for  the 


solution  X  is 


<Ad(X) 


ictqci 


— II  1  -  ,AdT.,Ad, 

k(E)k(ET)  d  d 


(4.27) 


where  Ad  -  A  -  B(R+BTXB)_i (ATXB+S)T. 

Proof:  With  reasoning  similar  to  that  for  the  continuous -time  problem, 
we  consider  a  perturbation  in  F,j(x)  on  the  order  of  a  perturbation  in 
the  constant  term,  i.e., 


isF.i  <  ic act  •  e 

a  —  in 


(4.28) 
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VV 


;^iv- 


and  the  desired  result  follows  from  Theorem  4.5. 

He  remark  here  that  when  E*I,  the  requirement  (4.25)  becomes  1  > 


lA^  IIA^I.  This  Is  always  true  for  "some"  norm  since  A^  is  the  closed- 

loop- system  matrix  and  all  of  its  eigenvalues  have  magnltutde  less  than 
unity.  Therefore,  its  spectral  radius  is  less  than  one  and  there 
always  exists  a  norm  that  is  arbitrarily  close  to  the  spectral  radius 
p(A^)  ([15],  Theorem  6.3.8).  Thus,  one  could  argue  in  this  case  that 


.  ictqci 

< A  (X)  -  - 21 - 

IXI(l-p2(AJ)) 


and  the  term  l-pz(A{j)  would  have  the  effect  of  raising  the  condition 

number  when  the  spectrum  of  the  closed -loop  system  has  a  pole  near  the 
unit  circle.  For  the  case  E  #  I,  (4.25)  can  be  a  very  restrictive 
requirement  and  is  a  subject  of  continuing  research. 

4.3  Balancing  to  Improve  Condition 

Considerable  success  has  been  attained  in  increasing  the  numerical 
accuracy  of  the  solutions  to  linear  equations  and  eigenvalue  problems 
by  appropriate  scaling  of  problem  parameters  [39],  [40].  Therefore,  it 
is  reasonable  to  expect  that  some  sort  of  balancing  or  problem  scaling 
could  improve  the  accuracy  of  the  numerical  computations  invblved  in 
computing  the  solution  to  the  GARE  and,  hence,  the  overall  “condition” 
of  the  problem.  One  could  argue  that  since  the  numerical  solution 
depends  fundamentally  on  the  QZ  algorithm  for  the  generalized  eigen- 
problem  [41],  that  balancing  the  eigensystem  for  the  generalized 
eigenproblem  XLy  ■  My  could  improve  the  Riccati  solution <  Ward  [42] 
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has  proposed  a  balancing  algorithm  specifically  designed  to  precede  QZ- 
type  algorithms.  This  balancing  consists  of  permutations  and  two  sided 
diagonal  transformations.  The  strategy  is  to  scale  L  and  M  so  chat 
their  elements  have  magnitudes  as  close  to  unity  as  possible.  Thus, 
the  large  elements  of  L  and  M  cannot  mask  the  effect  of  the  small 
elements,  as  can  often  be  the  case. 

Those  interested  in  the  details  of  the  algorithm  are  referred  to 
[42]  where  the  scaling  strategy  is  discussed.  Numerical  results  for 
example  elgenproblems  are  also  given  in  [42].  Numerical  experience  in 
solving  the  GARE  has  shown  that  this  balancing  strategy  can  increase 
the  accuracy  of  the  solution;  an  example  is  given  in  Chapter  5.  This 
balancing  can  also  increase  the  reliability  of  k(Vjj)  as  an  indicator 
of  condition  of  the  Riccatl  problem.  However,  one  does  pay  the  price 
in  Increased  computational  time,  about  10Z  of  the  QZ  algorithm 
computation  time. 

An  alternate  approach  to  the  direct  balancing  of  the  elgenproblem 
is  to  attempt  some  sort  of  coordinate  change  in  the  problem  which 
generates  the  Riccatl  equation.  That  is,  some  nonsingular  transforma¬ 
tion  T  to  change  coordinates  in  the  original  model  as  follows: 

x(t)  -  Tw(t)  (4.29) 

or 

^  -  Twfc  .  (4.30) 

There  is  reason  to  believe  that  the  coordinate  balancing  transformation 
proposed  by  Moore  [43]  may  result  in  a  transformed  Riccatl  equation 
that  is  more  amenable  to  numerics!  solution.  Moore  has  shown  [44]  that 


by  a  change  of  coordinates,  one  can  scale  the  relative  size  of  compon- 

T 

At  A*  t  T 

ents  of  e  o  and  e  C  .  Since,  as  we  have  already  pointed  out, 
stabillzablllty  of  the  model  influences  the  computational  accuracy  of 
the  Rlccatl  solution,  one  could  argue  that  a  change  of  coordinates 
increasing  the  relative  size  of  components  of  e^cB  would  be  benefi¬ 
cial  prior  to  computing  the  Riccati  solution.  Specifically,  Moore's 
’’internally  balanced”  coordinates  in  which  the  reachability  and  observ¬ 
ability  gramians  are  equal  and  diagonal  or  "input -normal”  coordinates 
where  the  reachability  gramian  is  identity  are  logical  choices. 

Laub  [6]  illustrated  the  effect  of  this  coordinate  type  balancing 
on  the  following  linear  optimal  control  problem:  Find  a  feedback 

controller  u(t)-Kx(t)  which  minimizes  the  performance  index 

J(u)  -  /"[xT(t)Ox(t)  +  uT(t)Ru(t)]dt 

with  plant  dynamics  given  by 

x(t)  -  Ax(t)  +  Bu(t) ;  x(0)  -  x0  . 

T  T 

Assume  Q»Q  >  0,  R-R  >0,  (A,B)  stabllizable  and  (A,C)  detectable,  where 
T 

C  C-Q  and  rank  (C)  *  rank  (Q).  Then  the  optimal  control  is  well  known 
to  be 

u(t)  -  -R"lBTXx 


where  X  solves  the  ARE 


Suppose  we  change  coordinates  via  (4.29).  Then  in  terms  of  the  new 
state  w(t),  our  problem  is  to  minimize 

/*[wT(t)TTQTw(t)  +  uT(t)Ru(t)]dt  (4.31) 

subject  to 

w(t)  -  (T-1 FT)w(t)  +  (T“lB)u(t).  (4.32) 

The  associated  solution  X,,  of  the  transformed  Riccati  equation  is 
related  to  the  original  X  by 

X  -  T~TX  T"1  (4.33) 

w 

One  can  see  from  (4.32)  and  (4.33)  that  if  T  is  ill-conditioned 
with  respect  to  inversion  and  the  balancing  is  applied  as  these  equa¬ 
tions  indicate,  then  this  technique  can  potentially  Introduce  more 
error  into  the  solution  than  would  originally  have  appeared.  Hence, 
the  opposite  of  the  intended  effect  could  occur.  Of  course,  one  could 
be  careful  to  only  choose  a  well- conditioned  T,  like  a  diagonal  scaling 
matrix  for  instance.  However,  if  T  is  a  modal  or  system  balancing  [43] 
transformation, it  could  be  quite  ill-conditioned. 

We  can  reduce  this  problem  significantly  in  our  generalized 
problem  formulation  framework  as  the  following  result  will  show. 

Theorem  4.7:  If  a  change  of  coordinates  of  the  form  (4.29)  is  made  in 
the  continuous -time  generalized  optimal  regulator  problem  of  section 
2.3.1  or  a  change  of  coordinates  of  the  form  (4.30)  is  made  in  the 
discrete-time  generalized  optimal  regulator  problem  of  Section  2.3.3, 
then  the  solution  X  to  the  original  problem  may  be  found  as  follows: 
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1)  replace  A  by  AT,  C  by  CT,  S  by  TT  S  and  E  by  ET  In  the 
appropriate  matrix  pencils  ((2.51)  or  (2.65)  for  the  continuous 
problem,  (2.67)  or  (2.68)  for  the  discrete  problem) 

2)  determine  the  appropriate  P,  Z  transformations  for  the 


modified  pencils  and  partition  Z  as  in  (2.61) 


3)  let  V 


El  0  - 
0  I  Z’ 


then  X  -  V2iVn~l 


(4.34) 


solves  the  original  GARE  ((2.30)  for  the  continuous  problem,  (2.48)  for 

T 

the  discrete  problem)  with  X  -  X  ^  0. 

Proof;  For  the  continuous-time  case,  the  problem  becomes  in  the  trans¬ 
formed  coordinates 


(4.35) 

(4.36) 


System:  ET  v(t)  »  ATW(t)  +  Bu(t) 

y(t)  -  CT  w(t) 

Criterion:  J  -  i-  /** (yTQy  +  uTRu  +  2wTTTSu)dt. 

1  0 

Application  of  Hamllton-Jacobi  theory  as  in  Section  2.3.1  gives  rise  to 
a  set  of  equations  to  which  we  can  associate  the  following  matrix 
pencil: 

(A-BR71 ST)T  -br-1bt 

I  -TT(CTQC-SR-lST)T  -T^A-BR^S1)1! 

This  pencil  can  also  be  obtained  by  performing  step  1  of  the  theorem  on 
the  pencil  associated  with  (2.57).  Now  we  can  factor  (4.37)  as 


ET 


T  T 

0  T  E 


(4.37) 


I  0 


0  T 


E  0 


0  E 


(A-BR_1ST) 


-br_1bt 


[_-(CTQC-SR-1ST)  -(A-BR-1ST)3j 


T  0 
0  I 


(4.38) 
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If  we  perform  the  ordered  transformation  on  (4.38)  by  finding  a  P  and  Z 
as  in  (2.60)  (as  stated  in  step  2  of  the  theorem)  and  compare  the 


result  with  (2.60),  we  find 


T  0  = 

0  1  Z  * 


Recall  from  Theorem  2.4  that 

x  -  v21vu-1  -  z21(ez11)-1 


(4.39) 


(4.40) 


Substitute  (4.39)  into  (4.40)  and  the  theorem  is  proved  for  the  contin¬ 
uous  case.  An  analogous  procedure  proves  the  discrete  case.  We  omit 
those  details. 

Note  from  the  theorem  that  the  numerous  occurrences  of  the  Inverse 
of  T  have  been  eliminated.  The  Inversion  of  T  is  essentially  required 
only  once  in  the  process  of  solving  the  linear  system 


xvn  -  v21 

for  the  Riccati  solution  X.  Computational  advantage  has  also  been 


realized  because  the  number  of  matrix  multiplications  required  in  the 
solution  process  has  been  reduced. 


•W 


CHAPTER  5 


NUMERICAL  EXPERIMENTS 

The  methods  for  solution  of  the  GARE  presented  in  Chapter  2,  the 
iterative  refinement  procedure  derived  in  Chapter  3,  and  the  condition 
estimates  for  the  computed  solution  of  Chapter  4  are  all  geared  to  the 
efficient  and  reliable  numerical  computation  of  the  Riccati  solution. 
This  chapter  describes  a  FORTRAN  software  package  (RICPACK)  developed 
to  aid  In  the  study  of  the  numerical  conditioning  of  the  algebraic 
Riccati  equation  and  other  closely  related  topics.  First,  a  brief 
description  of  the  package  will  be  given.  A  short  discussion  of  the 
algorithms  and  software  will  follow.  Finally,  numerical  examples  and 
results  will  be  given  to  illustrate  relevant  aspects  of  the  numerical 
solution. 

5.1  Software  Package  RICPACK 

RICPACK  was  developed  to  aid  in  the  study  of  the  numerical  condi¬ 
tioning  of  the  algebraic  Riccati  equation,  methods  for  improving  the 
condition,  and  iterative  refinement  of  the  solution.  The  FORTRAN 
subroutines  in  the  package  were  written  in  modular  form  and  are 
designed  to  facilitate  their  incorporation  in  some  larger  computer- 
aided  control  system  design  (CACSD)  package  or  computer-aided  systems 
and  control  analysis  and  design  environment  (CASCADE). 

A  FORTRAN  driver  program  has  also  been  written,,  primarily  as  a 
research  tool,  that  is  user-friendly  for  use  in  an  interactive  "term¬ 
inal”  type  environment.  The  program  prompts  for  all  necessary  input. 
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Convenient  input  default  options  exist  not  only  for  ease  of  data 
input,  but  also  for  exploitation  by  the  subroutines  to  reduce  the 
number  and  complexity  of  the  computations.  These  options  are  also 


designed  to  speed  the  input  of  more  “standard”  type  problems.  Table 
5.1  lists  the  default  options  available. 

TABLE  5.1 


Default  Options 


Input  Matrix 


Default  Value 


Identity 

Identity 


Identity  (or  enter 
diagonal  elements  only 
if  diagonal) 


Highlights  of  RICPACK  capabilities  include: 

(a)  Choice  of  calculation  of  the  stabilizing  (non-negative  defin¬ 
ite),  antl-stablllzing  (non-positive  definite),  or  just  any  (possibly) 
indefinite  or  nonsymmetrlc  solution  to  the  GARE. 

(b)  Coordinate  or  system  balancing  of  the  system  model  [43],  [45]; 
i.e.,  a  special  coordinate  transformation  on  (2.18)  or  (2.40)  of  the 


x(t)  -  T  w(t) 


■  T  w. 


(5.1) 


(5.2) 


respectively,  such  that  the  observability  and  reachability  gramians  of 
the  transformed  system  are  equal. and  diagonal. 


»  *  A  *  _!***■  ■  •  *  ■  *  *******  »  J  *  *  *  I  ’  k  «  •  *  »  *  * 


I  Residual I < 


•%  \ 
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(k)  Condition  estimates  for  the  Rlccatl  problem  derived  In  Chapter  4. 

5.2  Algorithms  and  Software 

Most  of  the  algorithmic  computations  (i.e.,  those  beyond  matrix 
algebra)  performed  In  RICPACK  employ  proven  stable  algorithms  coded 
Into  reliable,  portable  FORTRAN  software.  Two  sources  for  the  software 
are  LINPACK  [39]  and  EISPACK  [40],  [46].  Modified  LINPACK  software  Is 
used  for  linear  equation  solving  (the  BLAS  are  replaced  with  in-line 
code),  estimating  the  condition  of  a  non-diagonal  cost  matrix  R, 
solving  the  system  XVjj  »  Vji  and  In  estimating  the  condition  of  V^. 

The  LINPACK  singular  value  decomposition  (SVD)  is  used  In  compressing 
the  (2n-fm)  x  (2n+m)  pencil  when  necessary.  EISPACK-type  QZ  software  Is 
used  for  performing  the  QZ  transformation  and  generalized  eigenvalue 
calculations.  Other  EISPACK  software  is  used  in  calculating  the 
coordinate  balancing  [43],  [45]  transformation.  In  addition  to  software 
based  on  the  Bartels-Stewart  algorithm  [47]  for  Lyapunov  equations. 

The  Bartels-Stewart  based  software  is  also  used  in  condition  estimation 
and  Newton's  iteration  calculations.  Ward's  [42]  software  is  used  for 
the  eigenproblem  balancing,  and  Van  Dooren's  software  [44]  is  used  for 
ordering  the  transformed  pencil.  A  description  of  each  subroutine  used 
is  given  in  the  Appendix. 

5.3  Numerical  Examples 

The  following  simple  continuous -time  example  was  used  to  illus-  ' 
trate  the  numerical  properties  of  RICPACK  when  stabllizability  is  the 
key  factor: 
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Example  1 


1  ■  D  -°] x  +  [»] " 

y-  [1  I]x  (5.4) 

minimize  J^( y^y  +  uTu)dt  (5.5) 

This  system  is  stabilizable  for  e  *  0  and  completely  reconstructible. 
The  applicable  ARE  is 

ATX  +  XA  -  XBBTX  +  CTC  -  0  .  (5.6) 


The  “true"  solution  for  X  can  be  hand  calculated  for  comparison  pur¬ 
poses  as 


iwTt ^  i 

e2  2+/W2 


2+/I+?  *  4(2+/l+e2)2 


(5.7) 


Note  that  as  e-*-0  the  system  approaches  unstabllizabillty  and  the  (1,1) 
element  of  X  tends  to  Infinity. 

The  solution  to  this  problem  was  numerically  computed  on  a  DEC 
KL-10  (under  TOPS-20)  in  double  precision  using  RICPACK.  The  machine 
precision  is  near  10~18  in  this  case.  The  results  of  interest  are 
summarized  in  Table  5.2.  The  only  measure  of  condition  included  in  the 
table  la  the  condition  of  with  respect  to  inversion  because  the 

other  measures  did  not  give  an  indication  that  the  solution  accuracy 
was  degenerating  as  e+0.  We  note  here  that  the  data  in  Table  5.2  and 
the  succeeding  tables  that  are  expressed  as  a  power  of  10  are  rounded 
to  the  nearest  power  of  10. 
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♦Accuracy  In  correct  significant  digits. 


Some  useful  observations  can  be  made  on  this  data.  One  can  see 
that  for  this  example,  k(V^)  and  the  residual  (r)  are  both  good 
indicators  of  the  numerical  accuracy.  Since  machine  precision  is  near 
10”18,  one  would  expect  about  17  correct  digits  for  a  well-conditioned 
problem  (which  is  the  case  for  e*l).  The  data  indicates  that  one  digit 
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of  accuracy  is  lost  for  each  power  of  10  change  in  ic(Vn)  and  r.  This 
is  desirable  behavior  of  a  condition  estimate.  Note  that  Ward 
balancing  improves  the  condition  of  Vj  j  and  reduces  the  value  of  r  for 
the  same  value  of  e.  Ward  balancing  enables  solution  calculation  for 
smaller  values  of  e,  but  the  accuracy  is  not  as  smooth  a  function  of 
•c(Vn).  However,  the  residual  is  still  a  good  indicator  of  accuracy. 
Note  that  in  all  cases  with  a  reasonable  starting  guess  a  few  itera¬ 
tions  of  Newton ’s  method  restores  full  accuracy.  The  generalized 
eigenvalue  solution  was  used  as  a  starting  guess  and  was  considered 
reasonable  if  tc(V^)  <  1. /(machine  precision).  When  this  condition  was 
not  satisfied,  the  Newton  Iteration  failed  to  converge  to  the  desired 
solution. 

The  above  example  illustrates  that  stablllzablllty  of  the  model 
does  indeed  Influence  the  numerical  accuracy  of  the  Rlccati  solution. 
Also,  tc(Vji )  and  r  are  good  indicators  of  solution  accuracy  as  the 
model  approaches  unstabllizabllity.  However,  <(Vj ^ )  may  not  be  a  good 
Indicator  in  other  situations  as  the  following  example  will  show: 


Example  2 

-e  1  0  0 

1 

• 

-1  -e  0  0 

1 

X  ■ 

0  0  c  1 

x  + 

1 

0  0  -1  e 

Lu 

y  •  [l  1  1  l]x  (5.8) 

minimize  /"  (y*y  +  uTu)dt  (5.9) 

0 

The  model  is  completely  controllable  and  observable.  The  open- loop 
poles  are  at  ±e  ±j,  and  the  applicable  ARE  is 


1 
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A*X  +  XA  -  XBBTX  +  CTC  -  0  .  (5.10) 

The  solution  to  this  problem  was  numerically  computed  on  a  UNIVAC 
1100/83  In  double  precision  using  RICPACK.  The  machine  precision  Is 
near  10  In  this  case.  The  results  of  Interest  are  summarized  In 
Table  5.3.  Although  an  exact  hand  solution  was  not  possible  In  this 
case,  the  behavior  of  the  residual  can  be  used  to  judge  the  solution 
accuracy.  This  example  was  designed  to  assess  the  effect  of  the 
separation  of  the  closed -loop  spectrum  on  solution  accuracy  and  the 
ability  of  condition  estimates  to  detect  degrading  accuracy.  The 

column  CLP  In  Table  5.3  Indicates  the  position  (real  part)  of  the 
closed -loop  pole  nearest  the  Imaginary  axis  In  the  complex  plane. 


TABLE  5.3 

Numerical  Results  for  Example  2,  e  -  10~N 


N 

CLP 

<(vn) 

kAc(X) 

(4.18) 

•cB(X) 

(4.4) 

r 

(5.3) 

Newton 

Iterations 

r 

(5.3) 

0 

10° 

101 

10° 

102 

10-16 

1 

10-16 

3 

o 

1 

0* 

10° 

106 

107 

io-12 

1 

10“14 

5 

10-10 

10° 

1010 

1011 

10"8 

1 

io- 12 

7 

►— 

O 

1 

►— 

F 

10° 

1018 

1015 

io-8 

2 

10-16 

8 

10-16 

10° 

1016 

1017 

io-1 

2 

10-16 

9 

10-18 

10° 

1018 

1018 

IO- 1 

- 

- 

One  can  see  from  this  example  that  <(Vj})  provides  no  indication 
of  loss  of  accuracy  In  the  solution.  However,  in  this  example  kAc(X) 
and  ieB(X)  correlate  directly  with  the  behavior  of  the  residual,  and 
thus,  the  solution  accuracy.  Ward  balancing  of  the  eigenproblem  had  no 
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noticeable  effect  on  solution  accuracy  for  a  given  value  of  e.  Newton 
iterations  did  significantly  improve  solution  accuracy,  as  measured  by 
the  residual,  until  the  condition  KAg(X)  ■  <B(X)  -  1. /(machine 

precision).  At  this  point,  the  Iterations  failed  to  converge  to  the 
desired  solution,  as  was  the  case  in  example  1. 

This  example  shows  that  separation  of  the  closed -loop  spectrum 
does  Indeed  Influence  the  numerical  accuracy  of  the  Rlccati  solution 
and  that  the  condition  estimates  in  which  separation  is  a  factor 
provide  good  indicators  of  solution  degeneracy. 

The  following  example  illustrates  the  effect  of  111-conditioning 
of  the  R  weighting  matrix,  with  respect  to  Inversion,  on  the  numerical 
solution  for  the  continuous -time  case.  Recall  that  ill-conditioning  of 
R  is  not  necessarily  a  problem  in  the  discrete-time  case  since  its 
inverse  is  not  explicitly  required. 


atx  +  XA  -  XBR“1BTX  +  CTC  -  0  . 


(5.13) 
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The  solution  to  this  problem  was  numerically  computed  on  a  UNIVAC 
1100/83  in  double  precision  using  RICPACK.  The  machine  precision  is 
near  10"18  in  this  case.  The  results  of  interest  are  summarized  in 
Tables  5.4  and  5.5.  Ward  balancing  was  employed  in  the  calculations 
for  Table  5.4,  and  coordinate  balancing  was  employed  for  Table  5.5. 
Results  of  calculations  where  no  balancing  was  applied  were  nearly 
identical  to  those  of  Table  5.5  for  coordinate  balancing. 

TABLE  5.4 


Numerical  Results  for  Example  3,  Ward  Balancing,  e*10"N 


N 

•c(R) 

k(Vh) 

II 

1 

<B(X) 

(4.4) 

■  r 
(5.3) 

Newton 

iterations 

(5^3) 

n 

101 

101 

10° 

102 

>—* 

o 

1 

1 

10“17 

H 

102 

102 

102 

106 

\ 

o 

•—4 

1 

10“17 

D 

104 

103 

103 

© 

O 

H 

10-14 

2 

io-14 

6 

106 

104 

103 

10n 

10-12 

10 

10" 12 

8 

108 

105 

103 

1013 

10-10 

10 

10-1° 

10 

10l° 

106 

103 

1015 

10"8 

10 

00 

1 

o 

12 

1012 

107 

103 

1017 

10"6 

10 

10"8 

14 

1014 

108 

103 

1019 

10"5 

10 

10"5 

16 

1016 

109 

103 

10* 1 

10-3 

10 

10"3 

TABLE  5.5 


ftimerical  Results  for  Example  3.  System  Balancing,  e-10“N 


N 

<(R) 

«<VU) 

•cAc(X) 

(4.18) 

<B(X) 

(4.4) 

r 

(5.3) 

Newton 

Iterations 

(5%) 

0 

101 

103 

10° 

103 

10-15 

2 

IO" 17 

2 

102 

103 

102 

LO6 

10-15 

2 

io-17 

4 

104 

103 

103 

109 

io-12 

7 

10“14 

6 

106 

103 

103 

1011 

io-12 

10 

10-12 

8 

108 

103 

103 

1013 

10“9 

10 

10-11 

10 

IO10 

103 

103 

1015 

10-6 

10 

10-9 

12 

1012 

103 

103 

1017 

10“6 

10 

10-6 

14 

1014 

103 

103 

1019 

IO"2 

10 

H- • 

O 

1 

* 

16 

1016 

103 

103 

1021 

IO-1 

10 

10-2 

The  data  in  Table  5.4  indicate  that  k(R)  with  respect  to  inversion 
accurately  reflects  the  behavior  of  the  residual.  Also,  icCV^)  with 
respect  to  inversion  is  too  optimistic  in  its  estimation  of  the  problem 
condition  and  kB(X)  is  too  pessimistic.  kAc(X)  provides  no  informa¬ 
tion  in  this  case.  A  maximum  of  10  Newton  iterations  was  allowed,  and 
where  10  appears  in  the  table,  convergence  did  not  occur.  The  value 
for  the  residual  in  that  case  is  the  residual  associated  with  the 
solution  at  the  tenth  iteration.  One  can  see  that  the  ill -conditioning 
of  R  begins  to  dominate  the  numerical  accuracy  when  k(R}  >  106 .  Since 
R”1  is  Involved  in  the  Newton  iteration  calculations,  iterative 
improvement  does  not  improve  accuracy  when  <(R)  dominates.  The 


NWC  TP  6521 


convergence  criteria  used  for  the  Newton  Iterations  do  not  recognize 
this  fact  and  stop  the  Iterations. 

The  data  In  Table  5.5  are  essentially  the  same  as  that  in  Table  5.4 
except  for  one  important  aspect.  (c(Vtl)  with  respect  to  Inversion 
provides  no  information  on  the  problem  conditioning  once  k(R)  >  k(Vii). 
This  may  indicate  that  Ward  balancing  will  cause  other  sources  of 
problem  ill -conditioning,  beside  unstabillzability  of  the  model  and 
singularity  of  E,  to  be  reflected  in  kCV^). 

The  preceding  examples  illustrate  that  none  of  the  potential 
measures  of  conditioning  are  reliable  indicators  by  themselves. 
However,  numerical  experience  to  date  has  shown  that  in  all  cases  at 
least  one  of  the  measures  will  detect  the  degeneracy  of  numerical 
accuracy  as  it  occurs. 
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CHAPTER  6 


SECOND -ORDER  MODELS 


The  behavior  of  many  physical  systems  In  engineering  can  be 


modeled  by  the  following  system  of  equations: 
Mx(t)  +  (D  +  G)x(t)  +  Kx(t)  -  f(t) 


(6.1) 


where  x,f  e  Rn  and  M,D,G  and  K  e  R11*11.  Moreover,  in  describing 
physical  systems  one  can  make  the  following  assumptions  without  loss  of 


generality: 


M  -  Mt  >  0,  generalized  capacitive  storage 
D  ■  DT  >  0,  generalized  energy  dissipators 


(6.2) 


G  -  -CT, 


generalized  conservative  elements 


K  *  kT  0,  generalized  inductive  storage. 

The  model  (6.1)  can  describe  electrical,  mechanical,  thermal,  and  other 
systems  by  appropriate  choice  of  "through”  variables,  f (current  In 
electrical  systems,  force  in  mechanical  systems,  etc.)  and  "across” 
variables,  x(voltage,  displacement,  etc.)  [48] •  Analogies  exist  among 
the  various  types  of  systems. 

The  model  (6.1)  can  result  directly  from  lumped  parameter  models, 
or  finite  approximations  to  distributed  parameter  systems  described  by 
partial  differential  equations.  One  large  class  of  systems  of'  current 
importance  are  large  space  structures  (LSS),  which  are  large  distribu¬ 
ted  parameter  systems  that  are  most  often  discretized  by  the  finite 
element  method  into  the  form  (6.1)  [49].  The  problem  of  controlling 
LSS  motivated  the  studies  of  this  chapter,  although  the  results  are 
applicable  to  any  system  described  by  (6.1).  The  first  section  defines 
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the  LSS  framework  and  explores  the  inherent  structure  of  (6.1)  and 
(6.2)  that  might  be  computationally  exploitable.  The  second  section 
considers  specifically  the  solution  of  the  GARE  associated  with  second- 
order  models.  In  the  third  section,  criteria  are  discussed  for  the 
determination  of  controllability,  stabilizability ,  observability  or 
detectability  of  (6.1). 

6.1  Second-Order-Model  Structure  in  the  LSS  Framework 

In  the  LSS  framework,  x  is  a  displacement  vector,  f  is  a  force 
vector,  M  is  the  mass  matrix,  K  is  the  stiffness  matrix,  D  is  the 
damping  matrix,  and  the  G  matrix  gives  rise  to  gyroscopic  forces.  In 
general,  n  is  initially  of  very  high  order  and  the  aforementioned 
matrices  may  be  sparse.  The  force  vector  f  is  of  the  form 


f(t)  -  Fu(t) 


(6.3) 


where  u  e  Rm  is  the  control  input,  and  F  e  Rnxm  is  the  input  mat¬ 
rix.  One  usually  considers  an  output  equation  of  the  form 


where 


y(t)  -  Px(t)  +  V*(t) 


y  c  Rr  and  P,  V  e  Rrxn  • 


(6.4) 


The  traditional  method  for  dealing  with  models  of  the  form  (6.1) 
is  to  transform  to  the  equivalent  standard  matrix  first-order  (state 
variable)  form 


z(t)  **  Az(t)  +  Bu(t) 


(6.5) 


• . - , > 1 1 .  * 
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Other  first-order  realizations  of  this  form  which  are  of  potential 
Interest  Include: 


CH-G  M 
M  0 


D+G  M 
-M  0 


-K 

0 


[-] 


z  + 


-K  0 
-M 


z  + 


0  -K 
-K  -(D+G) 


] 


z  + 


(6.9) 

(6.10) 


(6.11) 


Note  that  (6.9)  and  (6.11)  have  symmetric  E  and  A  matrices  when  G  -  0, 
while  (6.1)  has  skew-symmetric  E  and  symmetric  A  whan  0*0.  Theme 
properties  are  computationally  advantageous  In  the  problems  In  the 
remainder  of  this  chapter. 

6.2  The  GARE  for  Second-Order  Models 


Consider  the  model  (6.7),  (6.4)  In  the  context  of  the  continuous- 


time  optimal  regulator  problem  of  section  2.3.1  and  let  D  -  D+G.  The 
GARE  for  this  problem  from  (2.30)  is 


-  D  s] 


i 


and  the  optimal  control  law  Is  given  by 


«°  ■  -irl[°  fTi  *  [5  «] 


-d*  \l  °1  *  \fl] 

q[p  v] 

-d3  Lo  mJ  LvtJ 

■’>  ■  [j  :] 

(6.12) 

e 

(6.13) 

employing  the  method  of 

section 

2.3,  then  one  is  faced  with  a  4n  x  4n  generalized  eigenproblem.  Clearly, 
this  is  undesirable  If  n  is  Itself  a  large  number. 
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Instead  we  expand  X  ■  X*  _>  0  Into  nxn  blocks  as  [50] 


*W  Y 

>  0,  W 

T 

Ya  Z 


WT  >  0,  WW+T 


Y,  Z  -  YTW+Y  >  0. 


(6.14) 


Substituting  (6.14)  into  (6.12)  and  oultiplylng  out  yields  the  follow¬ 


ing  three  equations 


0  -  -  YK  -  KYT  -  YFR-lFTYT  +  PTQP 
0  -  -  MZD  -  DTZM  -  MZFR-1FTZM  +  VTQV  +  HYT  +  YM 
0  -  W-YfJ-KZM  -  YFR-lFTZM  +  PTQV 


(6.15) 

(6.16) 
(6.17) 


We  see  that  (6.15)  is  an  nxn  Rlccatl-type  equation  from  which  Y  could 
be  obtained.  Once  Y  is  determined,  (6.16)  becomes  an  nxn  GARE  from 
which  we  can  solve  for  Z.  Equation  (6.17)  then  simply  defines  U. 

If  only  the  feedback  (6.13)  is  desired,  we  see  (by  substituting 
(6.14)  into  (6.13)  and  expanding)  that 


u°  -  -R"1FT(YTx  +  ZMx)  .  (6.18) 
Therefore,  one  need  not  solve  (6.17)  for  W. 

It  is  computationally  advantageous  to  only  require  the  solution  of 
two  nxn  Rlccatl  equations  rather  than  one  2n  x  2n  Rlccati  equation, 
since  the  solution  of  an  NxN  Rlccati  equation  takes  0(N3)  operations. 
To  realize  even  greater  computational  advantage,  we  consider  two  common 
special  cases.  First,  we  consider  the  case  when  there  is  no  damping  or 


gyroscopic  forces,  i.e.,  5*0.  One  can  see  that  (6.16)  is  then  no 
longer  a  GARE;  it  reduces  to  a  simpler  quadratic  equation. 
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Next,  consider  the  velocity  output  case,  l.e.,  P  -  0.  Then,  a 
solution  for  (6.15)  is  Y  ■  0.  One  can  solve  the  GARE  (6.16)  for  Z  and 
the  optlaal  control  becomes 

u°  -  -R-1FTZMx  (6.19) 

Obviously,  solving  one  nxn  GARE  is  computationally  simpler  than  solving 
a  2n  x  2n  GARE  regardless  of  the  size  of  n.  Based  on  these  results,  we 
can  state  the  following  theorem. 

Theorem  6.1:  Given  the  system 

M  x  (t)  +  D  x(t)  +  K  x(t)  -  F  u(t )  (6.20) 

with  output 

y(t)  -  Vx(t)  (6.21) 

and  criterion 

J  -  /  "(yT0y  +  uTRu)dt.  (6.22) 

0 

Then  the  unique  stabilizing  control  which  minimizes  (6.22)  is  given  by 

u°(t)  -  -R~1FTZMx(t)  (6.23) 

T 

where  Z  -  Z  ^  0  satisfies  the  GARE 

0  -  -  MZD  -  DTZM  -  MZFR~lFTZM  +  VTQV  .  (6.24) 

Proof:  The  theorem  is  a  restatement  of  the  results  of  this  section. 

6.3  Controllability  and  Observability  Criteria  for  Second- Order  Models 
We  now  seek  to  establish  controllability  and  observability  cri¬ 
teria  for  the  model  (6.1),  (6.3)  and  (6.4).  Controllability  and 

observability  of  this  model  have  been  shown  to  provide  important 
insights  into  modal  behavior  of  the  system  and  to  furnish  information 
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on  the  number  and  positioning  of  sensors  and  actuators  [ 51 ] ,  [52]. 

Also,  controllability  and  observability  information  can  be  used  In 
determining  which  modes  to  retain  when  performing  model  reduction  [53]. 

One  could  apply  the  traditional  methods  for  determining  control¬ 
lability  and  observability  to  the  transformed  first-order  model  (6.6) 
with  state  dimensions  2n.  Unfortunately,  the  resulting  standard  tests 
(cast  in  terms  of  a  2n-th  order  "A”  matrix)  do  not  take  advantage  of 
the  symmetry,  definiteness  and  sparsity  structure  of  the  matrices  H,  0, 
G,  and  K.  Also,  computational  problems  may  exist  if  H  is  near  singular 
or  n  is  very  large. 

Other  conditions  have  been  derived  by  Hughes  and  Skelton  [51], 
which  exploit  the  specialized  structure  of  (6.1)  and  (6.2)  for  the 
cases  D  •  0  and  DfG  ■  0.  However,  these  conditions  could  suffer  from 
computational  difficulties  since  they  require  knowledge  of  the  full 
modal  transformation  matrix,  whose  columns  are  the  eigenvectors  corre¬ 


sponding  to  the  eigenvalues  X  of  the  generalized  elgenproblem 
[XM  +  K]x  -  0 


(6.25) 


where  ±Xx  is  the  modal  frequencies.  An  equivalent  form  of  (6.25) 


is  the  simple  elgenproblem 


-M-1  Kx  -  Xx 


(6.26) 


since  M  is  nonsingular.  But  this  form  is  computationally  undesirable 
for  the  reasons  already  mentioned.  In  addition,  computation  of  the 
eigenvectors  of  (6.25)  and  (6.26)  is  ill-conditioned  whenever  the  X  is 
repeated  or  nearly  equal  [15],  which  is  often  the  case  in  LSS. 
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The  remainder  of  this  section  focuses  on  conditions  for  controlla¬ 
bility  (or,  more  generally,  stabillzability)  and  observability  (or, 
more  generally,  detectability)  which  take  advantage  of  the  structure  of 
(6.1)  and  (6.2),  but  extend  the  results  of  [ 51 ]  and  are  computationally 
more  tractable.  Most  of  the  computational  attractiveness  of  the  new 
criteria  accrue  from  the  fact  that  an  initial  modal  transformation  is 
not  necessary.  Thus,  if  just  a  few  "important**  modes  are  known — and 
there  exist  techniques  to  determine  just  selected  modes,  e.g.,  [54], 
[55] — these  modes  can  be  tested  for,  say  controllability,  by  a  test 
involving  just  the  model  matrices  M,  D,  G,  K, and  F. 

Definition  6.1: 

ft  :■  [Xi  :  X^  is  a  generalized  eigenvalue  of  the  problem 

r°  1 1-  x  [x  °1  ] 

l-K  -(D+G)j  |_0  MJ  ' 

:«  (modes  of  the  system  (6.7)} 

Also,  D+  :«  (X  e  a  :  Re  X  >  0} 

Q*  :»  (X  t  fl  :  Re  X  <  0) 

For  M  nonsingular  there  are  2n  modes,  and  ft+  and  ft~  are  the  sets  of 
the  unstable  and  stable  modes,  respectively.  Alternatively,  ft  can  be 
determined  in  terms  of  the  generalized  eigenvalue  problems 
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: 

4 


-K 


0 

-K  -(D+G) 


-  X 


-  X 


3- 


D+G  M 
M  0 


D+G 

-M 

-K 

0 


0 

M 


(6.30) 

(6.31) 

(6.32) 


(6.27) 

(6.28) 
(6.29) 


.1* 


1 
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whichever  yields  the  greatest  computational  advantage.  Computation  of 
eigenvalues  for  problems  of  the  form  (6.32)  Is  discussed  in  [54]  for 
the  case  G  -  0. 


Theorem  6.2:  (Hautus  [56])  The  system 


x  -  Ax  +  Bu  ;  x(t)  e  R 


(6.33) 


a)  Controllable  (stabilizable)  if  and  only  if 

rank  [A  -  XI,  B]  -  n;  for  all  X  e  A  (A)  (A+(A)) 

b)  Observable  (detectable)  if  and  only  if 


_(A  -  XI)  "  Q  ; 


for  all  X  e  A (A)  (A+(A)), 


where  A(A)  :■  (X:  (A-  XI)x  ■  0,  x*0] 


:■  Spectrum  of  A 


(6.34) 


(6.35) 


(6.36) 


Proof :  [56] 


Clearly  then,  the  system  (6.8)  with  E  nonsingular  is  controllable 
(stabilizable)  if  and  only  if 

rank  [A  -  XE,  B]  -  2n  for  all  X  e  A(E~lA)  (A+(E-LA))  (6.37) 


and  observable  (detectable)  if  and  only  if 


_A  -  XE  *  20 


for  all  X  e  A(E~1A)  (A+(E~lA)) 


(6.38) 


We  now  exploit  the  structure  of  A,  B,  C  to  derive  controllability, 
etc.,  criteria  directly  in  terms  of  M,  D,  G,  K,  F,  P,  and  V. 

Theorem  6.3:  The  system  (6.1),  ^6.3)  is  controllable  (stabilizable)  if 


and  only  if 


rank  [x2M  +  X(D+G)  +  K,  F]  -  n;  for  all  X  t  a  (fl+)  (6.39) 
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m 


;/v 

:;>J 


Proof :  By  Theorem  6.2, the  system  (6.1).  (6.3)  is  controllable  (stabil- 
izable)  if  and  only  if 

2n  ■  rank  [A  -  XE,  B] ;  for  all  X  e  ft  (G+) 


-  rank 


-  rank 


-  rank 


-XI  I  oj 
-K  -D-G-XM  Fj 

XM+D+G  ij  T -XI 

L  1  °JL-k 

[X2MfX(D-G)+K  0 
0  I 


from  (6.7)  and  (6.8) 


-XI  I 

-K  -D-G-XM 


ir  T 

0  -10  0 

F  -XI  I  0 

JL  0  0  Ij 


Clearly  this  obtains  if  and  only  if 

rank  [x2M+X(D+G)+R,  F]  -  n  for  all  X  e  n(n+). 

Note  that  X  is  a  scalar,  so  that  sparsity  in  the  problem  is  preserved. 
Also,  no  Inverses  and  no  Initial  transformations  are  necessary. 
Finally,  note  that  each  mode  of  the  system  can  be  checked  Individually 
without  transforming  the  system  to  modal  coordinates. 

Theorem  6.4:  The  system  (6.1),  (6.4)  is  observable  (detectable)  if  and 
only  if 


X2M+X(D+G)+K 


■  n  for  all  X  e  Jl  (0+) . 


(6.40) 


Proof;  By  Theorem  6.2  the  system  (6.1),  (6.4)  is  observable  (detect¬ 
able)  if  and  only  if 


“j  -*  ’ALa  >%j»  j' 


NWC  TP  6521 


2n  ■  rank 


■  rank 


■  rank 


C 

A  -  XE 


;  for  all  X  e  ft  (ft+) 


P  V 

-XI  I 

-K  -D-G-XM 


from  (6.7)  and  (6.8) 


rank 


-v  oir 

0  -XM-D-G  -I  1 1  - 

0  I  oj[_-l 

XV+P  0 j 

X2M+X(D+G)+K  0 

0  ij 


P  V 

XI  I 

K  -D-G-XM 


Clearly,  chl8  obtains  if  and  only  if 

XV+P 

X2M+X(D+G)+K 


rank 


-  n  for  all  X  e  0  (ft*). 


Note  that  an  alternative  proof  of  Theorems  6.3  and  6.4  is  to  observe 
that  the  results  are  essentially  restatements  of  the  spectral  criteria 
for  left  or  right  coprimeness  of  the  appropriate  polynomial  matrices 
(quadratic,  in  this  case).  However,  we  have  exploited  here  the  numer¬ 
ically  useful  characterization  of  ft  in  Definition  6.1,  and  the  proofs 
are  direct  and  require  no  polynomial  matrix  theory.  Several  special 
cases  of  Theorems  6.3  and  6.4  are  of  interest  in  many  systems  and  are 
now  stated  as  corollaries. 

Corollary  6.3.1:  When  D4G  •  0  (i.e.,  no  damping  or  gyroscopic  forces), 
the  system  (6.1),  (6.3)  is  controllable  if  and  only  if 

rank  [-w2M+K,  F]  -  n  ;  i  ■  1 . n  (6.41) 

where  -  /  X^  ;  X^  c  ft(M_1K)  (Note:  X^  >0). 

Proof :  When  D+G  ■  0 

ft  -  ;  i  -  l,...,n} 
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and  (6.41)  follows  directly  from  Theorem  6.3. 

Corollary  6.3.2:  When  D+G  -  0  and  the  system  (6.1),  (6.3)  is  in  modal 
form,  then  (6.1),  (6.3)  is  controllable  if  and  only  if 

rank  Fp  •  nf;  (r  -  1,...,R)  (6.42) 

where  Fr  are  partitioned  rows  of  the  modally  transformed  F  matrix 
corresponding  to  the  multiplicities  of  the  u>^;  n^  +•  ...  +  nR  - 
n,  and  R  is  the  number  of  distinct  o>£.  This  is  Theorem  1  of  [5l]. 

Proof :  In  modal  form 

M  -  I,  K  -  dlag[u>i2 , . .  .oi2  ],  F^  -  modally  transformed  F  mat¬ 
rix.  Then  (6.42)  follows  directly  from  Corollary  6.3.1. 

Corollary  6.3.3:  When  K  -  0  the  system  (6.1),  (6.3)  is  controllable  if 
and  only  if  rank  F  -  n. 

Proof :  When  K  *  0,  X-Oeft  and  the  corollary  follows  directly  from 

Theorem  6.3. 

Similarly,  we  can  state 

Corollary  6.4.1:  When  D+G  ■  0  the  system  (6.1),  (6.4)  is  observable  if 
and  only  if 

*  n  j  i  ■  1(, . . . , n .  (6. 43) 

Proof :  This  result  follows  directly  from  the  proof  of  Corollary  6.3.1 
and  Theorem  6.4. 

Corollary  6.4.2:  When  D+G  «  0  and  V  *  0  (i.e.,  no  rate  feedback)  the 
system  (6.1),  (6.4)  is  observable  if  and  only  if 
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rank 


P 

-w^MfK 


n;  i  -  1,.. . ,n. 


(6.44) 


Proof :  Follows  directly  from  Corollary  6.4.1. 

Corollary  6.4.3:  When  D+G  ■  0  and  the  system  (6.1),  (6.4)  is  in  modal 
form,  then  (6.1),  (6.4)  is  observable  if  and  only  if 

rankfjw^+P^]  -  n (r  -  1,...,  R)  (6.45) 

where  [ja>iVr  +  Pr]  are  the  suitably  partitioned  columns  of  the 
modally  transformed  [j(i>iV  +  P]  matrix. 

Proof :  In  modal  form  M  «  I,  K  -  diag[ui* ,...  ,<u*] ,  [jco^V  +  P]  -  modally 
transformed  [jwjV  +  p]  matrix  and  (6.45)  follows  directly  from 


(6.43). 


Corollary  6.4.4:  When  IHC  *  0,  the  system  (6.1),  (6.4)  is  in  modal 
form,  and  V  -  0  then  (6.1),  (6.4)  is  observable  if  and  only  if 
rank  P^  ■  n^;  (r  -  1,...,R).  (6.46) 

Proof:  Follows  directly  from  Corollary  6.4.3 

Corollary  6.4.5:  When  K  »  0  the  system  (6.1),  (6.4)  is  observable  if 
and  only  if  rank  P  -  n. 

Proof:  Set  w^  *  0  and  K  ■  0  in  (6.43). 

Note  that  the  above  Corollaries  have  obvious  analogues  in  termb  of 
stablllzabllity  and  detectability,  as  appropriate. 
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CHAPTER  7 
CONCLUSION 

We  have  examined  some  numerical  Issues  regarding  Che  solution  of  a 
very  general  form  of  the  algebraic  Riccati  equation  in  both  the  contin¬ 
uous-  and  discrete-time  formulations.  These  generalized  equations 
resulted  from  control  and  filtering  problems  for  systems  In  generalized 
state  space  form  with  performance  criteria  that  Included  cross-coupling 
between  the  state  and  input.  The  basic  solution  method  considered  was 
the  Schur  technique.  The  Schur  technique  was  preferred  because  of  Its 
good  numerical  properties,  especially  when  some  closed- loop-system 
eigenvalues  were  closely  spaced.  The  generalized  elgenproblem  frame¬ 
work  employed  allows  solutions  when  E  and  R  are  ill-conditioned  with 
respect  to  Inversion  without  undue  influence  of  the  ill -conditioning  on 
the  solution  process.  Singular  R  matrices  were  permissible  In  the 
discrete-time  case. 

A  Newton-type  Iterative  refinement  procedure  for  the  Riccati 
solution  was  derived.  In  the  most  general  case.  It  required  solution 
of  a  Sylvester  equation  at  each  Iteration.  Numerical  results  Indicated 
that  the  Iterative  refinement  improved  numerical  accuracy  significantly 
in  all  cases  where  accuracy  was  lost  due  to  ill-conditioning  of  the 
Riccati  problem  except  when  R  was  ill-conditioned  with  respect  to 
inversion  in  the  continuous-time  problem.  Lack  of  improvement  in  this 
case  was  attributed  to  the  fact  that  the  iterative  procedure  required 
R”1  explicitly  at  each  iteration. 


Condition  estimates  for  the  Riccati  problem  were  examined.  New 
condition  estimates  were  also  derived.  The  behavior  of  these  estimates 
and  their  ability  to  detect  degraded  accuracy  of  the  Riccati  solution 
were  evaluated  in  numerical  examples.  All  numerical  examples  employed 
the  software  package  RICPACK  for  the  solution  of  the  Riccati  equation. 
R1CPACK  was  developed  as  a  research  tool  to  aid  In  the  study  of  numeri¬ 
cal  Issues  related  to  the  solution  of  algebraic  Riccati  equations.  It 
was  found  that  no  single  condition  measure  considered  herein  would 
reliably  reflect  the  accuracy  of  the  Riccati  solution  in  all  cases. 
However,  ill-conditioning  was  always  detected  by  at  least  one  of  them. 

The  special  structure  of  linear  system  models  In  second-order  form 
was  considered.  Computational  advantage  was  sought  In  solving  Riccati 
equations  related  to  these  models.  Special  computational  advantage  was 
realized  in  the  velocity  feedback  control  problem.  Controllability  and 
observability  tests  for  second-order  models  were  derived  directly  in 
terms  of  the  system  model  matrices.  These  tests  had  the  additional 
advantage  of  enabling  one  to  test  Individual  model  modes . 

Some  topics  of  continuing  research  in  these  areas  will  Include: 

1.  Iterative  refinement  procedure  in  the  continuous-time  formula¬ 
tion  that  does  not  require  explicit  inversion  of  the  R  matrix. 

2.  Further  development  of  condition  estimates  for  the  Riccati 
problem,  especially  for  the  case  E  *  I  and  considering  other 
forms  for  the  perturbation  5F. 
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3.  Further  exploitation  of  the  special  structure  of  models  in 
second-order  form  In  all  types  of  control  and  filtering  compu¬ 
tations  . 

We  note  here  that  all  further  research  in  the  area  of  numerically 
solving  algebraic  Riccatl  equations  can  benefit  from  the  availability 
of  RICPACK. 
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The  purpose  of  this  appendix  is  to  provide  more  detailed  informa¬ 
tion  on  the  software  package  RICPACK  than  is  appropriate  for  the  main 
body  of  the  thesis.  However,  it  is  not  a  complete  documentation 
package  for  RICPACK.  First  there  will  be  a  discussion  of  the  hierarchy 
of  the  software  routines  employed  in  the  package.  Then  a  sample 
terminal  session  for  a  particular  example  problem  is  illustrated. 
Finally,  some  appropriate  software  listings  are  Included. 

Figure  A.  1  illustrates  the  hierarchy  of  the  software  routines. 
That  is,  the  routines  on  the  upper  levels  employ  the  routines  of  the 
lower  levels.  At  the  lowest  level  we  have  the  basic  matrix  manipula¬ 
tion  routines  like  add,  subtract,  multiply,  etc.,  and  some  simple 
combinations  of  these  basic  operations.  The  next  level  consists  of 
standard  routines  for  linear  equations,  eigenvalues  and  singular  value 
decomposition  (SVD).  Most  routines  in  this  level  are  from  LINPACK  or 
E1SPACK,  or  are  slight  modifications  to  routines  from  LINPACK  and 
EISPACK.  Subroutines  that  are  modified  have  the  modifications  noted  in 
the  comment  documentation  included  in  the  subroutine.  A  list  of  the 
Level  0  and  1  routines  is  given  in  Table  A.  1.  The  SVD  routine  is 
listed  separately  with  the  BLAS  routines  that  it  requires  from  LINPACK, 
as  it  is  the  only  routine  to  require  the  BLAS  and  could  be  modified  to 
eliminate  said  BLAS.  No  further  documentation  of  these  routines  is 
provided  in  this  appendix. 
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TABLE  A. 1 


I 

I 

t 

t  Subroutine  List  for  Levels  0  and  1 

♦ 


Level  0 

U 

*vel  1 

DSVDC 

BCORBK 

BALANC 

MLINEQ 

DAXPY 

D1NRM 

ORTHES 

DDOT 

DGECOM 

ORTRAN 

DNRM2 

muL 

DGEFAM 

QZHESW 

DROT 

DGESLM 

QZITW 

DROTG 

MQF 

QZVAL 

DSCAL 

MQFWO 

DSTSLV 

REBAKB 

DSWAP 

MSCALE 

ELMHES 

REDUCE 

MSUB 

GIV 

REDUC2 

MULA 

GRADBK 

ROTC 

MDLB 

GRADE Q 

ROTR 

MULWOA 

HQR 

SCALBK 

MJLWOB 

HQRORT 

SCALEG 

PERMUT 

EMTQL2 

SYMSLV 

SAVE 

LINEQ 

TRED2 

SEQOIV 

TRNATA 

TRNATB 

Subroutines  at  Levels  2  and  3  perform  more  specialized  tasks.  The 
task  title  is  given  in  Figure  A.l,  and  the  main  subroutine  performing 
the  task  is  shown  in  parenthesis.  The  major  comment  documentation  from 
these  main  subroutines,  as  it  appears  in  the  software,  is  given  in  this 
appendix.  This  documentation  can  be  consulted  for  more  details  on  the 
purpose  of  each  routine,  or  for  a  list  of  subroutines  that  each  of 
these  main  subroutines  may  call. 

Of  course,  the  main  driver  program  is  at  the  highest  level  (4),  and 
a  complete  listing  is  provided.  This  main  program  would  be  rewritten, 
or  at  least  extensively  modified  for  most  applications.  This  driver 
was  written  as  a  research  tool  and  as  such  performs  calculations  not 
relevant  to  many  analysis  and  design  applications.  A  higher  level 
language  would  be  more  appropriate  to  Interface  the  Level  0  through 
Level  3  subroutines  with  a  larger  CACSD  package  and  perform  the  neces¬ 
sary  input  and  output. 

A  sample  output  listing  from  an  actual  terminal  session  to  solve 
one  particular  problem  is  given  in  Figure  A. 2.  The  problem  is  Example 
I,  unbalanced  case,  for  e  -  .0001.  The  session  illustrates  most  of  the 
features  and  options  of  RICPACK.  The  computer  prompt  for  data  is  the 
”>”  sign.  The  listing  is  self-explanatory  and  will  not  be  elaborated 
further . 

The  main  program  listing  and  heading  documentation  comment  list¬ 
ings  for  the  main  subroutines  of  Levels  2  and  3  follow  Figure  A. 2. 
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@XQT  RICPACK. 

GENERALIZED  ALGEBRAIC  MATRIX  RICCATI  EQUATION  SOLVER 

ENTER  SYSTEM  ORDER  'N';  =  FOR  CONTINUOUS  TIME  PROBLEMS,  -  FOR  DISCRETE 
MAXIMUM  ORDER  =  20 

>2 

N  =  2 

ENTER  NUMBER  OF  SYSTEM  INPUTS  'M' 

AND  NUMBER  OF  SYSTEM  OUTPUTS  *L" 

>1.1 

M  =  1  L  =  1 


ENTER  FLAG  FOR  DESIRED  SOLUTION: 
-1  FOR  STABILIZING  SOLUTION 
0  FOR  ANY  SOLUTION 
♦  1  FOR  DESTABILIZING  SOLUTION 

>-1 

IORD  a  -1 


ENTER  BALANCING  FLAG:  0  FOR  WARD  BALANCING 

1  FOR  CO-ORDINATE  BALANCING 

2  FOR  NO  BALANCING 


>2 


IBAL  =  2 

NO  BALANCING  ATTEMPTED 


DO  YOU  WISH  TO  ITERATE  FOR  ROBUSTNESS  RECOVERY  (Y  OR  N) 
>N 

KFLAG  =  N 


DO  YOU  WISH  TO  ENTER  AN  *E*  MATRIX  (Y  OR  N) 
DEFAULT  IS  E  =  IDENTITY  MATRIX 


FIGURE  A. 2  Sample  terminal  session  output  listing. 
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>N 
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EFLAG  =  N 

USING  DEFAULT:  'E*  =  IDENTITY 

ENTER  THE  2X  2  SYSTEM  MATRIX 'A' BY  ROWS 

>1.0 

>0,-2 

THE  'A'  MATRIX  IS: 

1  00000000000000000  .000000000000000000 
000000000000000000  -2.00000000000000000 

ENTER  THE  2X  1  INPUT  MATRIX  "B'  BY  ROWS 

>  0001 
>0 

THE  *8"  MATRIX  IS: 

.  1 00000000000000000-003 
000000000000000000 

ENTER  THE  IX  2  OUTPUT  MATRIX  "C*  8Y  ROWS 

>1.1 

THE  "C  MATRIX  IS: 

1 .00000000000000000  1 .00000000000000000 

DO  YOU  WISH  TO  ENTER  A  CONTROL  WEIGHTING  MATRIX  ''  R"  (Y  OR  N) 

DEFAULT  IS  R  =  IDENTITY  MATRIX 

>N 

RFLAG  =  N 

USING  DEFAULT:  *R*  =  IDENTITY 

DO  YOU  WISH  TO  ENTER  AN  OUTPUT  WEIGHTING  MATRIX  *Q*  (Y  OR  N) 

>N 

FIGURE  A. 2  Sample  terminal  session  output  listing  (continued) 
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QFLAG  =  N 

USING  OEFAULT:  'Q'  =  IDENTITY 

DO  YOU  WISH  TO  ENTER  A  STATE/INPUT  CROSS-WEIGHTING  MATRIX  'S'  (Y  OR  N) 
DEFAULT  IS  S  =  ZERO  MATRIX 

>N 

SFLAG  =  N 

USING  DEFAULT:  'S'  =  ZERO  MATRIX 

THE  CLOSED  LOOP  EIGENVALUES  FOR  THE  STABILIZING  RICCATI  SOLUTION  ARE: 

(FOR  CONTINUOUS  TIME) 

-1 .00000000000000000  .000000000000000000 
-2.00000000000000000  .000000000000000000 

THE  STABILIZING  RICCATI  SOLUTION  IS: 

200000000.51 2931 183  .333333332798460259 

.33333333287 1560127  .249999999722222220 

X1N  =  .20000000  +  001  HIN  =  .20000000  +  001  K(Z1 1)  =  19402850  +  009 

KB(X)  =  .30000000  +  001  KA(X)  =  .50000000-008 

SEPEST  =  .20000000  +  001  SEP  =  .20000000  +  001 

DO  YOU  WISH  A  RESIDUAL  CALCULATION  (Y  OR  N) 

>Y 

RSFLG  =  Y 

RESIDUAL  1  -NORM/SOLUTION  1  -NORM  =  .129311 83873 1 697373-009 

DO  YOU  WANT  TO  SEE  THE  RESIDUAL  MATRIX  (Y  OR  N) 

>Y 

RMFLG  a  Y 


FIGURE  A. 2  Sample  terminal  session  output  listing  (continued) 
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THE  RESIDUAL  MATRIX,  BY  ROWS: 

•258623678237 1 99783-00 1  -.3205 1 5 1 6822S62 1316-010 
-.3205 1 5 1 6822562 1316-010  . 789603050697242 101-017 

DO  YOU  WANT  TO  TRY  ITERATIVE  IMPROVEMENT  (Y  OR  N) 

>Y 

NFLAG  =  Y 

ENTER  MAXIMUM  NUMBER  OF  NEWTON  ITERATIONS 
>5 

MAX.  NO.  OF  ITERATIONS  =  5 

CONVERGED  AFTER  2  ITERATIONS  OF  NEWTONS  METHOD 

CONVERGENCE  CRITERIA  1 

THE  REFINED  STABILIZING  SOLUTION  IS: 

200000000.499999999  .333333332777777778 

.333333332777777778  .249999999722222223 

THE  CLOSED  LOOP  EIGENVALUES  FOR  THE  STABILIZING  RICCATI  SOLUTION  ARE: 

•  1 .00000000499999998  .000000000000000000 

-2.00000000000000000  .000000000000000000 

RESIDUAL  1 -NORM/SOLUTION  1-NORM  *  .465661286234845640-017 

DO  YOU  WANT  TO  SEE  THE  RESIDUAL  MATRIX  <Y  OR  N) 

>Y 

RMFLG  =  Y 

THE  RESIDUAL  MATRIX,  BY  ROWS: 

.93 1 3225746 1 54785 1 6-009  .867361737988403547-0 1 8 
. 1 73472347597680709-017  -.2374427 1 4890844994-01 7 

THE  GAIN  MATRIX  IS: 

20000.0000499999998  .333333332777777779-004 

FIGURE  A. 2  Sample  terminal  session  output  listing  (continued). 
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C  THIS  IS  AH  INTERACTIVE  MAIN  DRIVER  PROGRAM  FOR  THE  SOFTWARE 

C  PACKAGE  RICPACK.  TIIIS  DRIVER  IS  FOR  SOLVING  THE  CONTINUOUS- 

C  TIME  OR  DISCRETE -TIME  GENERALIZED  OPTIMAL  REGULATOR  PROBLEM. 

C  THE  PROGRAM  PROMPTS  FOR  ALL  NECESSARY  INPUT.  FOR  THE  PROBLEM 
C  SPECIFICATIONS  AND  DETAILS  ON  THE  SOLUTION  METHOD,  SEE 
C 

C  REF.:  ARNOLD,  W.F.,  "ON  THE  NUMERICAL  SOLUTION  OF 
C  ALGEBRAIC  MATRIX  RICCATI  EQUATIONS,"  PHD  THESIS,  USC, 

C  DECEMBER  1983. 

C 

C  HISTORY: 

C  THIS  DRIVER  WAS  WRITTEN  BY  W.F.  ARNOLD,  NAVAL  WEAPONS  CENTER, 

C  CODE  35104,  CHINA  LAKE,  CA  93555,  AS  PART  OF  THE  SOFTWARE 

C  PACKAGE  RICPACK,  SEPTEMBER  1983. 

C 

C  SUBROUTINES  CALLED: 

C  BALANC,  BALCOR,  BCORBK ,  CMPRS,  ELMHES,  FBGAIN ,  HQR,  LYPCND , 

C  MADD,  MLINEQ ,  MOUT,  MQFWO,  MSCALE,  MSUB,  MULB,  NEWT,  RESID, 

C  RICSOL,  RINV,  SAVE,  SEPEST,  SEQUIV,  TRNATB 

C 

C  •••••FUNCTION  SUBPROGRAMS: 

DOUBLE  PRECISION  D1NRM 

INTEGER  IND(40) ,I,IBAL,INFO,  IORD,  J.  L,  LI,  L2.  L3,  M,  MAXIT, 
X  N,  NN,  NNPI ,  NNPJ ,  NNPM,  NOUT,  NPI,  NPJ,  NP1 ,  NR.  NRD,  NRT 
CHARACTER  CFLAG,  EFLAG,  KFLAG,  NFLAG,  QFLAG,  RDFLG ,  RFLAG, 

X  RMFLG,  RSFLG,  SFLAG 

DOUBLE  PRECISION  A(20,  20),  AS(20,  20),  B(20,  20),  C(20,  20), 

X  CQC ( 20 ,  20),  E(20,  20),  F<40,  40),  G(40,  40),  Q(20,  20), 

X  R(20 ,  20),  RS(20,  20),  S(20,  20),  U(60,  60),  WK(60,  20), 

X  Z(40,  40),  ALFI(40),  ALFR(40) ,  BETA ( 60 ) ,  CPERM(40), 

X  CSCALE ( 40 ) ,  CA,  CB,  CLT1N,  CLIN,  COND,  C1N ,  DG,  DGI ,  DGN, 

X  RSD,  RTOL,  SEP,  SR,  TEMP,  X1N 

LOGICAL  TYPE 

DATA  NR,  NRD,  NRT  /  20,  40,  60  / 

RTOL  s  1.0D  +  06 
INFO  =  0 
NOUT  =  6 
TYPE  =  .TRUE. 

CFLAG  =  'Y' 

WRITE  (NOUT,  800)  NR 
C 

C  READ  IN  PROBLEM  PARAMETERS 
C 

READ  (5,  •)  N 
WRITE  (NOUT,  802)  N 
IF  (N.GT.O)  GOTO  10 
TYPE  =  .FALSE. 

N  s  -  N 
10  CONTINUE 

IF  (N.LE.NR  .AND.  N.NE.O)  GOTO  20 


112 


NWC  TP  6521 


WRITE  (NOUT.  804) 

STOP 

20  CONTINUE 

WRITE  (NOUT,  806) 

READ  (5,  *)  M,  L 

WRITE  (NOUT,  808)  M,  L 

IF  (M.EQ.O)  WRITE  (NOUT,  810) 

IF  (L.NE.O)  GOTO  30 
WRITE  (NOUT,  812) 

STOP 

30  CONTINUE 

WRITE  (NOUT,  814) 

READ  (5,  •)  IORD 
WRITE  (NOUT,  816)  IORD 
IF  (M.EQ.O)  IORD  =  -  1 
WRITE  (NOUT,  818) 

READ  (5,  •)  IBAL 
WRITE  (NOUT,  820)  IBAL 

IF  (IBAL.NE.O  .AND.  IBAL.NE. 1 )  WRITE  (NOUT,  822) 
WRITE  (NOUT,  824) 

READ  (5,  834)  KFLAG 
WRITE  (NOUT,  826)  KFLAG 
IF  (KFLAG. NE.'Y')  GOTO  40 
WRITE  (NOUT,  828) 

READ  (5,  •)  DGI 
WRITE  (NOUT,  830)  DGI 
40  CONTINUE 
NN  =  2  *  N 
NNPM  =  NN  -*■  M 
IF  (IBAL.EQ. 1 )  GOTO  50 

READ  IN  PROBLEM  MODEL  MATRICES 

WRITE  (NOUT,  832) 

READ  (5,  834)  EFLAG 
WRITE  (NOUT,  836)  EFLAG 
IF  (EFLAG. EQ.'Y')  GOTO  60 
50  CONTINUE 

WRITE  (NOUT,  838) 

GOTO  80 
60  CONTINUE 

WRITE  (NOUT,  840)  N,  N 
DO  70  I  =  1,  N 

READ  (5,  *)  (E(I,  J),  J  =  1,  N) 

70  CONTINUE 

WRITE  (NOUT,  842) 

CALL  MOUT (NOUT,  NR,  N,  N,  E) 

80  CONTINUE 

WRITE  (NOUT,  844)  N,  N 
DO  90  I  =  1,  N 


READ  (5,  *)  (A(I,  J),  J  s  1,  N) 
90  CONTINUE 

CALL  SAVE (NR,  NR,  N,  N,  A,  AS) 

WRITE  (NOUT,  846) 

CALL  MOUT(NOUT,  NR,  N,  N,  A) 

IF  (M.EQ.O)  GOTO  110 
WRITE  (NOUT,  848)  N,  M 
DO  100  I  =  1,  N 

READ  (5,  »)  (B(I,  J).  J  =  1,  M) 
100  CONTINUE 

WRITE  (NOUT,  850) 

CALL  MOUT(NOUT,  NR,  N,  M,  B) 

110  CONTINUE 

WRITE  (NOUT,  852)  L,  N 
DO  120  I  =  1,  L 

READ  (5,  »)  (C(I,  J),  J  s  1,  N) 
120  CONTINUE 

WRITE  (NOUT,  854) 

CALL  MOUT (NOUT,  NR,  L,  N,  C) 

CALL  SAVE (NR,  NR,  L,  N,  C,  CQC) 

IF  (M.EQ.O)  GOTO  180 
WRITE  (NOUT,  856) 

READ  (5.  834)  RFLAG 
WRITE  (NOUT,  858)  RFLAG 
IF  (RFLAG. EQ.'Y*)  GOTO  130 
WRITE  (NOUT,  860) 

GOTO  180 
130  CONTINUE 

WRITE  (NOUT,  862) 

READ  (5,  834)  RDFLG 
WRITE  (NOUT,  864)  RDFLG 
IF  (RDFLG. NE. ’Y')  GOTO  160 
WRITE  (NOUT,  866)  M 
DO  150  J  =  1,  M 

DO  140  I  =  1,  M 

R(I,  J)  =  O.ODO 
140  CONTINUE 
150  CONTINUE 

READ  (5,  *)  (R(I,  I),  I  =  1,  M) 
WRITE  (NOUT,  868) 

CALL  MOUT (NOUT,  NR,  M,  M,  R) 

GOTO  180 
160  CONTINUE 

WRITE  (NOUT,  870)  M,  M 
DO  170  I  s  1,  M 

READ  (5,  •)  (R(I,  J),  J  =  1,  M) 
170  CONTINUE 

WRITE  (NOUT,  868) 

CALL  MOUT (NOUT,  NR,  M,  M,  R) 

180  CONTINUE 
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IF  (KFLAG.EQ. * Y * )  GOTO  230 
WRITE  (NOUT,  872) 

READ  (5,  834)  QFLAG 
WRITE  (NOUT,  874)  QFLAG 
IF  (QFLAG. EQ. *Y')  GOTO  190 
WRITE  (NOUT,  876) 

GOTO  210 
190  CONTINUE 

WRITE  (NOUT,  878)  L,  L 
DO  200  I  =  1,  L 

READ  (5,  *>  (Q(I,  J),  J  =  1,  L) 

200  CONTINUE 

WRITE  (NOUT,  880) 

CALL  MOUT(NOUT,  NR,  L,  L,  Q) 

210  CONTINUE 

FORM  THE  MATRIX  PRODUCT  CT*Q»C  AND  STORE  IN  CQC 

CALL  TRNATB(NR,  NRD,  L,  N.  C,  G) 

IF  (QFLAG. NE.'Y')  GOTO  220 

CALL  MULB(NR ,  NR,  L,  L,  N,  Q,  CQC,  CPERM) 

220  CONTINUE 

CALL  MULB(NRD,  NR,  N,  L,  N,  G,  CQC.  CPERM) 

GOTO  320 

FORM  ROBUSTNESS  RECOVERY  TERM 


230  CONTINUE 
ITER  =  1 

WRITE  (NOUT,  882) 

READ  (5,  834)  QFLAG 
WRITE  (NOUT,  874)  QFLAG 
IF  (QFLAG. EQ.’Y1)  GOTO  240 
WRITE  (NOUT,  884) 

GOTO  260 
240  CONTINUE 

WRITE  (NOUT.  886)  N,  N 
DO  250  I  s  1,  N 

READ  (5,  •)  (Q(I,  J),  J  s  1,  N) 

250  CONTINUE 

WRITE  (NOUT,  880) 

CALL  MOUT (NOUT,  NR,  N,  N,  Q) 

260  CONTINUE 

CALL  TRNATB(NR,  NRD,  L,  N,  C,  G) 

CALL  MULB(NRD,  NR,  N,  L,  N,  G,  CQC,  CPERM) 
CALL  SAVE (NR,  NRT,  N,  N,  CQC,  WK) 

CALL  MSCALE(NR,  N,  N,  DGI,  CQC) 

DG  =  DGI 

IF  (QFLAG. EQ.'Y')  GOTO  280 
DO  270  I  =  1 .  N 
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cqc(i,  i)  =  cqc(i,  d  ♦  i.odo 
270  CONTINUE 
GOTO  29C 
280  CONTINUE 

CALL  MADD(NR,  NR,  NR,  N,  N,  CQC,  Q,  CQC) 

290  CONTINUE 

DO  310  I  =  1,  N 

DO  300  J  =  1,  N 

WK( J+N,  I)  =  WK(J,  I) 

300  CONTINUE 
310  CONTINUE 
320  CONTINUE 

IF  (M.EQ.O)  GOTO  350 
WRITE  (NOUT,  888) 

READ  (5.  834)  SFLAG 
WRITE  (NOUT,  890)  SFLAG 
IF  (SFLAG. EQ. 'Y')  GOTO  330 
WRITE  (NOUT.  892) 

GOTO  350 
330  CONTINUE 

WRITE  (NOUT,  894)  N,  M 
DO  340  I  =  1,  N 

READ  (5.  •)  (S(I,  J),  J  s  1,  M) 

340  CONTINUE 

WRITE  (NOUT,  896) 

CALL  MOUT(NOUT,  NR,  N,  M,  S) 

350  CONTINUE 

CALCULATE  CO-ORDINATE  BALANCING  TRANSFORMATION  IF  REQUESTED 
IF  (IBAL.NE.1)  GOTO  360 

CALL  BALCOR(NR,  NRD,  L,  M,  N,  A,  B,  C,  CQC,  S,  E,  Q,  Z,  ALFI, 
X  ALFR,  BETA,  IND,  INFO,  SFLAG,  TYPE) 

IF  (INFO.EQ.O)  GOTO  360 
WRITE  (NOUT,  898) 

IF  (INFO.EQ. 1 )  WRITE  (NOUT,  900) 

IF  (INF0.EQ.-1)  WRITE  (NOUT,  902) 

IF  (INFO.EQ. -2)  WRITE  (NOUT,  904) 

WRITE  (NOUT,  906) 

READ  (5,  834)  RSFLG 
IF  (RSFLG. NE. 'Y')  STOP 
IBAL  =  0 
360  CONTINUE 

TAKE  INVERSE  OF  R  MATRIX  AND  STORE  IN  Q 

IF  (RFLAG.NE. 'Y')  GOTO  420 
IF  (RDFLG.NE.’Y')  GOTO  380 
DO  370  I  r  1,  M 

Q(I,  I)  s  I.ODO  /  R(I ,  I) 
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370  CONTINUE 
GOTO  420 
380  CONTINUE 

CALL  SAVE(NR,  NRD,  M,  M,  R,  G) 

DO  400  J  =  1,  M 

DO  390  I  s  1,  M 

Q(I,  J)  =  0.0D0 
390  CONTINUE 

Q(J,  J)  =  1.0D0 
400  CONTINUE 

CALL  MLINEQ(NRD,  NR,  M,  M,  G,  Q,  COND,  IND,  CPERM) 

WRITE  (NOUT,  908)  COND 
IF  (COND.LT.RTOL)  GOTO  410 
WRITE  (NOUT,  910) 

READ  (5.  834)  CFLAG 
IF  ( CFLAG. NE.’Y')  GOTO  430 
410  CONTINUE 

WRITE  (NOUT,  912) 

420  CONTINUE 

SET  UP  MATRIX  PENCIL  FOR  CONTINUOUS-TIME  PROBLEM  USING 
R  INVERSE 

IF  (SFLAG.EQ.'Y')  CALL 

X  SEQUIV (NR,  NRD,  M,  N,  A,  B,  CQC,  S,  Q,  F,  G,  RDFLG,  RFLAG ) 
CALL  RINV(NR,  NRD,  N,  NN,  M,  E,  A,  B,  CQC,  Q,  G,  F,  RS.  CPERM, 

X  RDFLG,  RFLAG,  EFLAG,  IBAL,  TYPE) 

GOTO  440 
430  CONTINUE 

SET  UP  MATRIX  PENCIL  FOR  CONTINUOUS-TIME  PROBLEM  USING 
COMPRESSION  TECHNIQUE  WHEN  R  IS  ILL-CONDITIONED  WITH 
RESPECT  TO  INVERSION 

WRITE  (NOUT,  914) 

IF  (.NOT. TYPE)  CFLAG  =  'Y* 

KFLAG  =  'N' 

CALL  CMPRS(NR,  NRD,  NRT,  N,  NN,  NNPM,  M,  E,  A,  B,  CQC,  R,  S,  G, 
X  F,  U,  WK,  CPERM,  CSCALE,  BETA,  EFLAG,  SFLAG,  IBAL,  INFO) 

IF  (INFO.NE.O)  GOTO  720 
440  CONTINUE 

IF  (TYPE)  GOTO  470 

CONVERT  PENCIL  TO  DISCRETE-TIME  PENCIL  IF  DISCRETE  PROBLEM 

NP1  s  N  ♦  1 
DO  460  J  =  NP1 ,  NN 
DO  450  I  s  1,  NN 
TEMP  =  G(I,  J) 

G(I,  J)  s  F(I ,  J) 
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F( I ,  J)  =  -  TEMP 
450  CONTINUE 
460  CONTINUE 
470  CONTINUE 

IF  (IBAL.EQ.1)  CALL  BCORBKCNR,  NRD,  M,  N,  A,  B,  CQC,  S,  Z) 

SAVE  PENCIL  IF  ITERATING  FOR  ROBUSTNESS  RECOVERY 

IF  (KFLAG.NE. 'Y')  GOTO  500 
CALL  SAVE (NRD,  NRT,  NN,  NN,  G.  U) 

DO  490  I  =  1,  N 
NPI  a  N  +  I 
NNPI  =  NN  +  I 
DO  480  J  =  1,  N 
NPJ  =  N  +  J 
NNPJ  =  NN  ♦  J 
U(NNPJ,  I)  s  F( J ,  I) 

U(NNPJ ,  NPI)  =  F( J ,  NPI) 

U(J,  NNPI)  =  F(NPJ ,  I) 

U (NPJ ,  NNPI)  =  F(NPJ ,  NPI) 

480  CONTINUE 
490  CONTINUE 
500  CONTINUE 

COMPUTE  THE  RICCATI  SOLUTION 

CALL  RICSOL(NR,  NRD,  NN,  N,  G,  F,  E,  Z,  ALFR ,  ALFI,  BETA, 

X  CPERM ,  CSCALE,  IND,  IORD,  IBAL,  TYPE,  EFLAG) 

IF  (IND(I).NE.O)  GOTO  690 
INFO  =  IND (2) 

IF  (CPERM(I).EQ.I.OD  +  20)  CALL  SAVE (NRD,  NRD,  N,  N,  Z,  F) 
IF  (IBAL.NE. 1 )  GOTO  530 

ESTIMATE  CONDITION  OF  CO-ORDINATE  BALANCING  TRANSFORMATION 
WITH  RESPECT  TO  INVERSION 

DO  520  J  *  1,  N 

DO  510  I  =  1,  N 

G(I,  J)  =  O.ODO 
510  CONTINUE 

G(J,  J)  =  1.0D0 
520  CONTINUE 

CALL  SAVE (NR,  NRD,  N,  N,  E,  Z) 

CALL  MLINEQ(NRD,  NRD,  N,  N,  Z,  G,  COND,  IND,  CSCALE) 

WRITE  (NOUT,  916)  COND 
530  CONTINUE 

OUTPUT  THE  SOLUTION 


WRITE  (NOUT,  918) 
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IF  (.NOT. TYPE)  GOTO  560 
WRITE  (NOUT,  920) 

DO  550  I  s  1,  NN 

IF  (ALFR(I) .GE.0.0D0)  GOTO  550 
IF  (BETA(I) .NE.O.ODO)  GOTO  540 
WRITE  (NOUT,  922) 

GOTO  550 
540  CONTINUE 

ALFR(I)  =  ALFR(I)  /  BETA(I) 

ALFI(I)  =  ALFI(I)  /  BETA(I) 

WRITE  (NOUT.  »)  ALFR(I),  RLFI(I) 

550  CONTINUE 
GOTO  580 
560  CONTINUE 

WRITE  (NOUT,  924) 

DO  570  I  *  1,  NN 

IF  (BETA(I) .EQ.O.ODO)  GOTO  570 
IF  (DABS(ALFRd))  .GE.BETA(I))  GOTO  570 
ALFR(I)  =  ALFR(I)  /  BETA ( I ) 

ALFI(I)  =  ALFI(I)  /  BETA(I) 

WRITE  (NOUT,  •)  ALFR(I) ,  ALFI(I) 

570  CONTINUE 
580  CONTINUE 

IF  (INFO.NE. 0)  GOTO  700 
IF  (CPERM(I).EQ.I.OD  +  20)  GOTO  710 
IF  (M.EQ.O)  WRITE  (NOUT,  926) 

IF  (I0RD.EQ.-1)  WRITE  (NOUT,  928) 

IF  (IORD.EQ.O)  WRITE  (NOUT,  930) 

IF  (IORD.EQ. 1 )  WRITE  (NOUT,  932) 

CALL  MQUT(NOUT,  NR,  N,  N,  F) 

C 

C  COMPUTE  CONDITION  ESTIMATES 
C 

X1N  =  D1NRM(NRD,  N,  N,  F) 

CALL  SAVE (NR,  NR,  N,  N,  CQC,  C) 

IF  (SFLAG.EQ.'Y')  CALL 

X  SEQUIV (NR ,  NRD,  M,  N,  AS,  B,  C,  S,  Q,  Z,  G,  RDFLG ,  RFLAG) 
C1N  =  D1NRM(NR,  H,  N,  C) 

CALL  FBGAIN (NR ,  NRD,  NRD,  N,  M',  A,  B,  E,  R,  Q,  S,  F,  Z,  G, 

X  CSCALE,  IND,  EFLAG,  RDFLG,  RFLAG,  SFUG,  TYPE) 

WRITE  (NOUT,  934)  X1N,  C1N,  CPERM(I) 

CALL  MULB(NR,  NRD,  N,  M,  N,  B,  Z,  CSCALE) 

CALL  MSUB(NR,  NRD,  NRD,  N,  N,  A,  Z,  Z) 

IF  (TYPE)  GOTO  600 

CLIN  =  D1NRM(NRD,  N,  N,  Z) 

CALL  TRNATA(NRD,  N,  Z) 

CLT1N  a  D1NRM(NRD,  N,  N,  Z) 

CALL  3ALANC(NRD,  N,  Z,  LI,  L2 ,  CSCALE) 

CALL  ELMHES(NRD,  N,  LI,  L2,  Z,  IND) 

CALL  HQR(NRD,  N,  LI,  L2,  Z,  ALFR,  ALFI,  L3) 


119 


.s'..-. 


■ 


IF  (L3.NE.0)  WRITE  (NOUT,  936) 

SR  =  0.0D0 
DO  590  I  =  1.  H 

SR  =  DMAXKSR,  DSQRT ( ALFR ( I ) * ALFR ( I )  +  ALFI(I)«ALFI(I)) ) 
590  CONTINUE 

TEMP  =  1.0D0  -  CLIN  •  CLT1N 
IF  (TEMP.LE. O.ODO)  TEMP  =  1.0D0  -  SR  •  SR 
CA  =  DFLOAT(N)  •  (C1N  /  X1N  +  2.0D0  *  CLIN  «  CLT1N)  /  TEMP 
WRITE  (NOUT,  938)  CA,  SR,  CLIN,  CLT1N 
GOTO  620 
600  CONTINUE 
LI  =  0 
L2  =  1 

CALL  LYPCND(NRD,  NRD,  N,  Z,  F,  G,  ALFR,  ALFI,  BETA,  LI,  L2) 
CALL  SEPESTCNRD,  N,  Z,  G,  SEP,  LI) 

IF  (L1.NE.0)  WRITE  (NOUT,  940) 

TEMP  =  DABS (ALFR ( 1 )) 

DO  610  I  r  2,  N 

TEMP  s  DMIN1 (TEMP,  DABS ( ALFR ( I ) ) ) 

610  CONTINUE 

CB  =  (C1N  /  X1N  +  2.0D0  •  D1NRM(NR,  N,  N,  AS)  +  X1N  • 

X  D1NRM(NR,  N.  N,  RS))  /  SEP 
CA  =  C1N  /  (X1N  •  SEP) 

TEMP  =  2.0D0  •  TEMP 
WRITE  (NOUT,  942)  CB,  CA,  TEMP,  SEP 
620  CONTINUE 

SET  UP  FOR  ANOTHER  ROBUSTNESS  RECOVERY  ITERATION 

IF  (KFLAG.NE. ’ Y 1 )  GOTO  680 
WRITE  (NOUT,  944) 

READ  (5,  834)  KFLAG 
WRITE  (NOUT,  826)  KFLAG 
IF  (KFLAG. NE.'Y’)  GOTO  650 
WRITE  (NOUT,  946) 

READ  (5,  •)  DGN 
WRITE  (NOUT,  830)  XN 
IF  (IBAL.EQ.1  .AND.  ITER. EQ. 1 ).  CALL 
X  MQFWO(NRT,  NR,  N,  WK,  E,  CPERM) 

ITER  =  ITER  ♦  1 
X  =  X  -  XN 
X  640  I  =  1,  N 
NPI  s  N  +  I 
NNPI  =  HN  +  I 
X  630  J  =  1,  N 
NPJ  =  N  +  J 
NNPJ  =  NN  ♦  J 

U(NPJ ,  I)  =  U(NPJ ,  I)  ♦  X  •  WK(J ,  I) 

F( J ,  I)  =  U(NNPJ ,  I) 

F( J  ,  NPI)  =  U (NNPJ ,■  NPI) 
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F(NPJ ,  I)  =  U(J,  NNPI ) 

F(NPJ ,  NPI)  =  U (NPJ ,  NNPI) 

630  CONTINUE 

640  CONTINUE 

CALL  SAV£(NRT,  NRD,  NN,  NN.  U.  G) 

DG  =  DGN 
GOTO  500 
650  CONTINUE 

DG  =  DG  -  DGI 
DO  670  I  =  1,  N 

DO  660  J  =  1,  N 

CQC(J,  I)  =  CQC( J ,  I)  ♦  DG  •  WK(J+N,  I) 
660  CONTINUE 

670  CONTINUE 
680  CONTINUE 

IF  ( CFLAG .  NE . '  Y  * )  STOP 


COMPUTE  THE  RESIDUAL 


WRITE  (NOUT,  948) 

READ  (5.  834)  RSFLG 
WRITE  (NOUT,  950)  RSFLG 
IF  (RSFLG. NE.’Y1)  GOTO  730 

CALL  RESID(NR,  NR,  NRD,  NRD,  N,  M,  E,  A,  B.  CQC,  R,  S,  Q,  RS 
X  F,  G,  Z.  CPERM,  IND,  RTOL,  EFLAG,  RFLAG,  RDFLG,  SFLAG, 

X  RSD,  TYPE,  NOUT) 

WRITE  (NOUT.  952)  RSD 
WRITE  (NOUT,  954) 

READ  (5,  834)  RMFLG 
WRITE  (NOUT,  956)  RMFLG 
IF  (RMFLG. NE.'Y*)  GOTO  730 
WRITE  (NOUT,  958) 

CALL  MOUT(NOUT,  NR,  N,  N,  RS) 

GOTO  730 

690  CONTINUE 

WRITE  (NOUT,  960) 

STOP 

700  CONTINUE 

WRITE  (NOUT,  962) 

STOP 

710  CONTINUE 

WRITE  (NOUT,  964) 

CALL  MOUT (NOUT ,  NRD,  N,  N,  F) 

STOP 

720  CONTINUE 

WRITE  (NOUT,  966) 

STOP 

730  CONTINUE 

IF  (I0RD.NE.-1)  STOP 


n  o  o  on 
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APPLY  NEWTON’S  METHOD  FOR  ITERATIVE  IMPROVEMENT 

WRITE  (NOUT,  968) 

READ  (5,  834)  NFLAG 
WRITE  (NOUT,  970)  NFLAG 
IF  (NFLAG. NE.'Y')  GOTO  750 
WRITE  (NOUT.  972) 

READ  (5,  •)  MAXIT 

WRITE  (NOUT,  974)  MAXIT 

CALL  SAVE(NRD,  NRT,  N,  N,  F,  U) 

CALL  NEWT(NR,  NR,  NRD,  NRT,  N,  M,  E,  A,  B,  CQC,  R,  S,  Q,  RS, 

X  U,  G,  F,  Z,  ALFR,  ALFI,  BETA,  IND,  RTOL,  EFLAG,  RFLAG , 

X  RDFLG,  SFLAG,  SEP,  TYPE,  MAXIT,  INFO,  NOUT) 

IF  (INFO.NE. 0)  GOTO  760 
WRITE  (NOUT,  976)  MAXIT 
WRITE  (NOUT,  978)  IND(1) 

CALL  RESID(NR,  NR,  NRD,  NRT,  N,  M,  E,  A,  B,  CQC,  R,  S,  Q,  RS, 
X  U,  G,  Z,  CPERM,  IND,  RTOL,  EFLAG,  RFLAG,  RDFLG,  SFLAG, 

X  RSD,  TYPE,  NOUT) 

WRITE  (NOUT,  980) 

CALL  MOUT (NOUT,  NRT,  N,  N,  U) 

WRITE  (NOUT,  918) 

DO  740  I  =  1,  N 

WRITE  (NOUT,  *)  ALFR (I),  ALFI (I) 

740  CONTINUE 

WRITE  (NOUT.  952)  RSD 

CALL  SAVE (NRT,  NRD,  N,  N,  U,  F) 

WRITE  (NOUT,  954) 

READ  (5.  834)  RMFLG 
WRITE  (NOUT,  956)  RMFLG 
IF  (RMFLG. NE.'Y')  GOTO  750 
WRITE  (NOUT,  958) 

CALL  MOUT (NOUT,  NR,  N,  N,  RS) 

750  CONTINUE 

CALCULATE  OPTIMAL  FEEDBACK  GAIN  MATRIX 

CALL  FBGAIN (NR ,  NRD,  NRD,  N,  M,  A,  B,  E,  R,  Q,  S,  F,  Z,  G, 

X  CPERM,  IND,  EFLAG,  RDFLG,  RFLAG,  SFLAG,  TYPE) 

IF  ( CPERM ( 1 ).GT. RTOL)  WRITE  (NOUT,  984)  CPERM(I) 

WRITE  (NOUT,  982) 

CALL  MOUT (NOUT,  NRD,  M,  N,  Z) 

STOP 

760  CONTINUE 

IF  (INF0.EQ.-3 )  GOTO  770 
WRITE  (NOUT,  986)  INFO 
CALL  MOUT (NOUT,  NRT,  N,  N,  U) 

WRITE  (NOUT,  988)  BETA( 1 ) 

CALL  RESID(NR,  NR,  NRD,  NRT,  N,  M,  E,  A,  B,  CQC,  R,  S,  Q,  RS, 
X  U,  G,  Z,  CPERM,  IND,  RTOL,  EFLAG,  RFLAG,  RDFLG,  SFLAG, 
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X  RSD,  TYPE,  NOUT) 

GOTO  740 
770  CONTINUE 

WRITE  (NOUT,  990) 

STOP 

800  FORMAT (/'  GENERALIZED  ALGEBRAIC  MATRIX  RICCATI  EQUATION  SOLVER' 
X  //,'  ENTER  SYSTEM  ORDER  "N":  +  FOR  CONTINUOUS  TIME  PROBLE ' 

X  'MS,  -  FOR  DISCRETE',  /,  '  MAXIMUM  ORDER  =',  13.  /) 

802  FORMAT  (/'  N  =' ,  13/) 

804  FORMAT  (//'  ORDER  EXCEEDS  MAXIMUM  ***  EXECUTION  TERMINATED',  /) 
806  FORMAT  (//'  ENTER  NUMBER  OF  SYSTEM  INPUTS  "M"',  /, 

X  '  AND  NUMBER  OF  SYSTEM  OUTPUTS  "L"'/) 

808  FORMAT  (/'  M  =' ,  13,  5X,  'L  =',  13/) 

810  FORMAT  (/'  SYSTEM  HAS  NO  INPUTS,  RICCATI  EQUATION  DEGENERATES', 
X  '  TO  A  LYAPUNOV  EQUATION'/) 

812  FORMAT  (//’  SYSTEM  HAS  NO  OUTPUT  •»»  EXECUTION  TERMINATED'/) 

814  FORMAT  (/'  ENTER  FLAG  FOR  DESIRED  SOLUTION:',  /, 

X  '  -1  FOR  STABILIZING  SOLUTION',/,'  0  FOR  ANY  SOLUTION',  /, 

X  '  +1  FOR  DESTABILIZING  SOLUTION'/) 

816  FORMAT  (/'  IORD  =  ',  12/) 

818  FORMAT  (/'  ENTER  BALANCING  FLAG:  0  FOR  WARD  BALANCING',  /, 

X  '  1  FOR  CO-ORDINATE  BALANCING'/, 

X  ’  2  FOR  NO  BALANCING'/) 

820  FORMAT  (/'  IBAL  =' .  12/) 

822  FORMAT  (/'  NO  BALANCING  ATTEMPTED'/) 

824  FORMAT 

X  (/'DO  YOU  WISH  TO  ITERATE  FOR  ROBUSTNESS  RECOVERY', 

X  '  (Y  OR  N)'/) 

826  FORMAT  (/'  KFLAG  =  ',  A1/) 

828  FORMAT  (/'  ENTER  INITIAL  GAMMA  PARAMETER'/) 

330  FORMAT  (/'  DG  =’ ,  D26.18/) 

832  FORMAT  (/’  DO  YOU  WISH  TO  ENTER  AN  "E”  MATRIX  (Y  OR  N)' ,  /, 

X  '  DEFAULT  IS  E  =  IDENTITY  MATRIX'/) 

834  FORMAT  (1A1) 

836  FORMAT  (/'  EFLAG  =  ',  A1/) 

838  FORMAT  (/'  USING  DEFAULT:  "E"  =  IDENTITY'/) 

840  FORMAT  (/'  ENTER  THE  ',  13,  '  X  ',  13,  ’  MATRIX  "E"  BY  ROWS'/) 
842  FORMAT  (/'  THE  "E"  MATRIX  IS:'/) 

844  FORMAT 

X  (/’  ENTER  THE  M3,'  X  M3,'  SYSTEM  MATRIX  "A"  BY  RpWS'/) 
846  FORMAT  (/'  THE  "A"  MATRIX  IS:'/) 

848  FORMAT 

X  (/'  ENTER  THE  ',13,'  X  M3,'  INPUT  MATRIX  "B”  BY  ROWS'/) 
850  FORMAT  (/'  THE  "B"  MATRIX  IS:'/) 

852  FORMAT 

X  (/'  ENTER  THE  ',  13,'  X  M3,'  OUTPUT  MATRIX  ”0"  BY  ROWS'/) 
854  FORMAT  (/'  THE  "C"  MATRIX  IS:'/) 

856  FORMATS ’  DO  YOU  WISH  TO  ENTER  A  CONTROL  WEIGHTING  MATRIX  "R"’ , 
X  '  (Y  OR  N)',  /,  '  DEFAULT  IS  R  =  IDENTITY  MATRIX'/) 

858  FORMAT  (/'  RFLAG  =  ',  A1/) 
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860  FORMAT  (/'  USING  DEFAULT:  "R"  a  IDENTITY'/) 

862  FORMAT (/'  DO  YOU  WISH  TO  ENTER  A  DIAGONAL  "R"  MATRIX (Y  OR  N)'/) 
864  FORMAT  (/•  RDFLG  a  A1/) 

866  FORMAT  (/'  ENTER  THE  13.  '  DIAGONAL  ELEMENTS  AS  A  ROW'/) 

868  FORMAT  (/'  THE  "R"  MATRIX  IS:'/) 

870  FORMAT  (/'  ENTER  THE  13.  '  X  '  ,13, '  INPUT  WEIGHTING  MATRIX’, 
X  *  "R"  BY  ROWS'/) 

872  FORMAT (/'  DO  YOU  WISH  TO  ENTER  AN  OUTPUT  WEIGHTING  MATRIX  •’Q"', 
X  *  (Y  OR  N)',  /,  '  DEFAULT  IS  Q  a  IDENTITY  MATRIX'/) 

874  FORMAT  (/'  QFLAG  a  ',  A1/) 

876  FORMAT  (/'  USING  DEFAULT:  "Q*  a  IDENTITY’/) 

878  FORMAT  (/'  ENTER  THE  ',  13,  *  X  ',  13, 

X  '  OUTPUT  WEIGHTING  MATRIX  "Q"’ ,  •  BY  ROWS'/) 

880  FORMAT  (/'  THE  "Q"  MATRIX  IS:'/) 

882  FORMAT  (/’  DO  YOU  WISH  TO  ENTER  A  STATE  WEIGHTING  MATRIX  "Q"' , 

X  '  (Y  OR  N)',  /,  •  DEFAULT  IS  Q  a  IDENTITY  MATRIX'/) 

884  FORMAT  (/'  USING  DEFAULT:  "Q"  a  IDENTITY'/) 

886  FORMAT  (/'  ENTER  THE',  13,  '  X  *,  13, 

X  '  STATE  WEIGHTING  MATRIX  "Q"' ,  '  BY  ROWS'/) 

888  FORMAT  (/'  DO  YOU  WISH  TO  ENTER  A  STATE/INPUT  CROSS-WEIGHTING ' , 
X  '  MATRIX  "S"  (Y  OR  N)' ,  /,  •  DEFAULT  IS  S  a  ZERO  MATRIX'/) 
890  FORMAT  (/'  SFLAG  a  ',  A1/) 

892  FORMAT  (/'  USING  DEFAULT:  "S"  a  ZERO  MATRIX'/) 

894  FORMAT  (/'  ENTER  THE  ',  13.  *  X  ',  13, 

X  '  CROSS -WEIGHTING  MATRIX  "S'",  '  BY  ROWS'/) 

896  FORMAT  (/'  THE  "S"  MATRIX  IS:'/) 

898  FORMAT  (/’  CO-ORDINATE  BALANCING  UNSUCCESSFUL  BECAUSE:') 

900  FORMAT  ('  CAN  NOT  COMPUTE  SOLUTION  TO  LYAPUNOV  EQUATION  /) 

902  FORMAT  ('  CONTROLLABILITY  GRAMMIAN  NOT  NUMERICALLY  ?D.'/> 

904  FORMAT  ('  CONVERGENCE  FAILURE  IN  EIGENVALUE  ROUTINE'/) 

906  FORMAT (/'  DO  YOU  WISH  TO  CONTINUE  WITH  WARD  BALANCING (Y  OR  N)'/) 
908  FORMAT 

X  (/'  INVERSION  CONDITION  ESTIMATE  FOR  "IP  MATRIX  a’,D26.18,/) 
910  FORMAT  (/'  DO  YOU  WISH  TO  CONTINUE  USING  "R"  INVERSE  (Y  OR  N) '/) 
912  FORMAT  (/'  PROCEEDING  WITH  SOLUTION  USING  "R"  INVERSE'/) 

914  FORMAT  (/'  PROCEEDING  WITH  COMPRESSION  TECHNIQUE'/) 

916  FORMAT  (/’  A  CONDITION  ESTIMATE  FOR  THE  CO-ORDINATE  BALANCING’, 

X  '  TRANSFORMATION  IS:',  /,  D26.18/) 

918  FORMAT  (/ 

X  '  THE  CLOSED  LOOP  EIGENVALUES  FOR  THE  STABILIZING  RlfCATI', 

X  ’  SOLUTION  ARE:’/) 

920  FORMAT  ('  (FOR  CONTINUOUS  TIME)'/) 

922  FORMAT  (/•  INFINITE  EIGENVALUE'/) 

924  FORMAT  ('  (FOR  DISCRETE  TIME)’/) 

926  FORMAT  (/’  THE  LYAPUNOV  SOLUTION  IS:’/) 

928  FORMAT  (/'  THE  STABILIZING  RICCATI  SOLUTION  IS:'/) 

930  FORMAT  (/'  A  RICCATI  SOLUTION  IS:'/) 

932  FORMAT  (/•  THE  DESTABILIZING  RICCATI  SOLUTION  IS'/) 

934  FORMAT 

X  (/'  X1N  a',D15.8,5X,'C1N  a',D15.8,5X,  ’K(ZII)  a’,  D15.8/) 
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936  FORMAT  (/'ERROR  IN  HQR  '/) 

938  FORMAT  (/'  KA(X)  ='.  D15.8,  10X,  'SR  =  ' ,  D15.8,  /,  '  CLIN  =  ', 

X  D15.8,  10X,  ' CLT1N  =• ,  D15.8/) 

940  FORMAT  (/'  ERROR  IN  SEPEST'/) 

942  FORMAT  (/’  KB(X)  =  ',  D15 .8, 10X, 'KA(X)  *',D15.8,/,'  SEPEST  =  ' , 

X  D15.8,  10X,  'SEP  =' ,  D15.8,  /) 

944  FORMAT  (/'  DO  YOU  WISH  TO  ENTER  ANOTHER  GAMMA  (Y  OR  N)'/) 

946  FORMAT  (/'  ENTER  NEW  VALUE  FOR  GAMMA'/) 

948  FORMAT  (/'  DO  YOU  WISH  A  RESIDUAL  CALCULATION  (Y  OR  N)'/) 

950  FORMAT  (/’  RSFLG  =  ',  A1/) 

952  FORMAT  (/'  RESIDUAL  1 -NORM/SOLUTION  1-NORM  =’ ,  D26.18,  /) 

954  FORMAT  (/'  DO  YOU  WANT  TO  SEE  THE  RESIDUAL  MATRIX  (Y  OR  N)'/) 
956  FORMAT  (/'  RMFLG  =  ',  A1/) 

958  FORMAT  (/'  THE  RESIDUAL  MATRIX,  BY  ROWS:'/) 

960  FORMAT  (/'  MORE  THAT  50  ITERATIONS  REQUIRED  BY  QZITW’ ,  /, 

X  '  «*  EXECUTION  TERMINATED'/) 

962  FORMAT  (/'  CONVERGENCE  FAILURE  IN  ORDER  ROUTINE',  /, 

X  '  »**  EXECUTION  TERMINATED'/) 

964  FORMAT (/'  SCHUR  VECTOR  MATRIX  SINGULAR  TO  WORKING  PRECISION',/, 
X  '  THEREFORE,  SOLUTION  BY  THIS  TECHNIQUE  NOT  POSSIBLE',  /, 

X  '  THE  SCHUR  VECTOR  MATRIX  IS:'/) 

966  FORMAT  (/'  COMPRESSION  FAILURE  EXECUTION  TERMINATED'/) 

968  FORMAT  (/'  DO  YOU  WANT  TO  TRY  ITERATIVE  IMPROVEMENT  (Y  OR  N)'/) 
970  FORMAT  (/'  NFLAG  =  ',  A1/) 

972  FORMAT  (/'  ENTER  MAXIMUM  NUMBER  OF  NEWTON  ITERATIONS'/) 

974  FORMAT  (/'  MAX.  NO.  OF  ITERATIONS  =’ ,  14/) 

976  FORMAT 

X  (/'  CONVERGED  AFTER  M3,'  ITERATIONS  OF  NEWTONS  METHOD'/) 
978  FORMAT  ('  CONVERGENCE  CRITERIA',  12/) 

980  FORMAT  (/'  THE  REFINED  STABILIZING  SOLUTION  IS:'/) 

982  FORMAT  (/'  THE  GAIN  MATRIX  IS:'/) 

984  FORMAT  (//'  R  +  GT»X«G  ILL-CONDITIONED  WRT  INVERSION’,  /, 

X  '  CONDITION  NUMBER  =' ,  D26.18,  //) 

986  FORMAT  (/'  NEWTON  ITERATION  FAILED,  INFO  r  ',  13/ , 

X  '  THE  SOLUTION  AT  THE  LAST  ITERATION  IS:’/) 

988  FORMAT  (/'  THE  1-NORM  OF  THE  ERROR  IS:',  D26.18/) 

990  FORMAT (/'  SORRY,  PROGRAM  NOT  ABLE  TO  PERFORM  NEWTON  ITERATION', 
X  /,  '  FOR  ARBITRARY  E -MATRIX  AT  THIS  TIME!!!'/) 

END 
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SUBROUTINE  BALCOR  (NR,NRZ,L,M,N,A,B,C,CQC,S,T,WK,Z,WK1 ,WK2, 
X  WK3 , IPVT , ISTAB , SFLAG , TYPE ) 


•••••PARAMETERS: 

INTEGER  NR , NRZ , L , M , N , IPVT ( N ) , ISTAB 
CHARACTER  SFLAG 

DOUBLE  PRECISION  A(NR,N) ,B(NR,M) ,C(NR,N) ,CQC(NR,N) ,S(NR,M) , 
X  T(NR,N),WK(NR,N),Z(NRZ,N),WK1 (N ) ,WK2 (N ) ,WK3 (N ) 

LOGICAL  TYPE 


•••••LOCAL  VARIABLES: 
INTEGER  I, J,NPI,NPJ 


•••••FORTRAN  FUNCTIONS: 
NONE. 


•••••SUBROUTINES  CALLED: 

BCORBK,  BLCRDC ,  BLCRDD,  MMUL,  MQFWO,  MULB,  TRNATB 


•••••PURPOSE: 

GIVEN  THE  MODEL  MATRICES  A,  B  AND  C  FOR  A  FIRST  ORDER  LINEAR 
MODEL  IN  STATE  SPACE  FORM,  THIS  SUBROUTINE  CALCULATES  A 
BALANCING  TRANSFORMATION,  T,  SUCH  THAT  IF  T  WAS  APPLIED  TO  THE 
MODEL  AS  A  CHANGE  OF  COORDINATES,  THE  REACHABILITY  AND 
OBSERVABILITY  GRAMM IANS  WOULD  BE  EQUAL  AND  DIAGONAL.  HOWEVER, 
T  IS  ACTUALLY  APPLIED  TO  THE  MODEL  MATRICES  IN  A  SPECIAL  WAY 
FOR  SOLUTION  OF  THE  OPTIMAL  REGULATOR  PROBLEM  BY  RICPACK. 


REF.:  ARNOLD,  W.F.,  "ON  THE  NUMERICAL  SOLUTION  OF 
ALGEBRAIC  MATRIX  RICCATI  EQUATIONS,"  PHD  THESIS,  USC, 
DECEMBER  1983. 


‘PARAMETER  DESCRIPTION: 


ON  INPUT: 


INTEGER 

ROW  DIMENSION  OF  THE  ARRAYS  CONTAINING  THE  A,  B, 
C,  CQC,  S,  T  AND  WK  MATRICES  AS  DECLARED  IN  THE 
MAIN  CALLING  PROGRAM  DIMENSION  STATEMENT; 


INTEGER 

ROW  DIMENSION  OF  THE  ARRAY  CONTAINING  THE  Z  MATRIX 
AS  DECLARED  IN  THE  MAIN  CALLING  PROGRAM  DIMENSION 
STATEMENT ; 


INTEGER 
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C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ROW  DIMENSION  OF  THE  C  MATRIX; 


M  INTEGER 

COLUMN  DIMENSION  OF  THE  B  AND  S  MATRICES; 


N  INTEGER 

ORDER  OF  THE  SQUARE  MATRICES  A,  CQC  AND  T 
COLUMN  DIMENSION  OF  THE  C  MATRIX; 


A  REAL (NR, N) 

MODEL  SYSTEM  MATRIX,  ALTERED  BY  THIS  ROUTINE; 


B  REAL(NR.M) 

MODEL  INPUT  MATRIX; 


C  REAL(NR.N) 

MODEL  OUTPUT  MATRIX,  ALTERED  BY  THIS  ROUTINE; 


CQC  REAL (NR, N) 

MATRIX  PRODUCT  CT*Q*C  WHERE  T  DENOTES  MATRIX 
TRANSPOSE,  ALTERED  BY  THIS  ROUTINE; 


S  REAL (NR, M) 

STATE  -  INPUT  CROSS-WEIGHTING  MATRIX.  ALTERED  BY 
THIS  ROUTINE; 

WK  REAL (NR, N) 

SCRATCH  ARRAY  OF  SIZE  AT  LEAST  N  BY  N; 


WK1 ,WK2 ,WK3 

REAL(N) 

WORKING  VECTORS  OF  SIZE  AT  LEAST  N; 


IPVT  INTEGER (N) 

WORKING  VECTOR  OF  SIZE  AT  LEAST  N; 


SFLAG  CHARACTER 

FLAG  SET  TO  'Y'  IF  S. IS  OTHER  THAN  THE  ZERO  MATRIX 


TYPE  LOGICAL 

=  .TRUE.  FOR  A  CONTINUOUS-TIME  MODEL 
=  .FALSE.  FOR  A  DISCRETE-TIME  MODEL. 


ON  OUTPUT: 

A  CONTAINS  THE  MATRIX  PRODUCT  F»T; 

CQC  CONTAINS  THE  MATRIX  PRODUCT  (C*T)-TRANSPOSE»Q»C*T; 
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S  CONTAINS  THE  MATRIX  PRODUCT  T-TRANSPOSE*S; 

T  REAL (NR, N) 

CONTAINS  THE  BALANCING  TRANSFORMATION,  T; 

Z  REAL(NRZ.N) 

2N  BY  (N+M)  MATRIX  WITH  THE  ORIGINAL  A  MATRIX  STORED 

IN  THE  UPPER  LEFT  N  BY  N  BLOCK,  CQC  IN  THE  LOWER  LEFT 
N  BY  N  BLOCK,  S  IN  THE  UPPER  RIGHT  N  BY  M  BLOCK  AND 
B  IN  THE  LOWER  RIGHT  N  BY  M  BLOCK; 

ISTAB  INTEGER 

ERROR  FLAG  WITH  THE  FOLLOWING  MEANINGS 
=  ZERO,  NORMAL  RETURN 

=  NONZERO,  A  BALANCING  TRANSFORMATION  COULD  NOT  BE 
CALCULATED  AND  THE  A,  CQC  AND  S  MATRICES  ARE 
RETURNED  UNALTERED. 

***** ALGORITHM  NOTES: 

NONE. 

•••••HISTORY: 

THIS  SUBROUTINE  WAS  WRITTEN  BY  W.F.  ARNOLD,  NAVAL  WEAPONS 
CENTER,  CODE  35104,  CHINA  LAKE,  CA  93555,  AS  PART  OF  THE 
SOFTWARE  PACKAGE  RICPACK,  SEPTEMBER  1983. 
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SUBROUTINE  RICSOL (NR , NRD , NN , N , G ,F , E , Z , ALFR , ALFI , BETA , CPERM , 

X  CSCALE , IND , IORD , IBAL , TYPE , EFLAG ) 

•••••PARAMETERS: 

INTEGER  NR , NRD , NN , N , IND (NN ) , IORD , IBAL 
CHARACTER  EFLAG 

DOUBLE  PRECISION  G(NRD,NN ) ,F (NRD.NN ) , E (NR ,N) ,Z(NRD,NN) , 

X  ALFR (NN ) ,ALFI (NN ) .BETA (NN ) .CPERM (NN ) .CSCALE (NN ) 

LOGICAL  TYPE 

•••••LOCAL  VARIABLES: 

INTEGER  I, IERR, IF AIL , IGH . J , LOW . NPJ 
DOUBLE  PRECISION  COND.EPS, EPS  1 
LOGICAL  MATZ 

•••••FORTRAN  FUNCTIONS: 

NONE. 

•••••SUBROUTINES  CALLED: 

BALGBK .  BALGEN.  MLINEQ,  MULB,  ORDER.  QZHESW.  QZITW,  QZVAL 

•••••PURPOSE: 

GIVEN  THE  2N  BY  2N  MATRIX  PENCIL 
LAMBDA  •  F  -  G 

THIS  SUBROUTINE  FINDS  THE  ORTHOGONAL  MATRIX  Z  SUCH  THAT 
Q  •  (LAMBDA  •  F  -  G)  •  Z 

IS  IN  GENERALIZED  ORDERED  REAL  SCHUR  FORM.  FURTHERMORE,  THE 
UPPER  LEFT  N  BY  N  BLOCK  OF  THE  TRANSFORMED  PENCIL  CONTAINS 
THE  EIGENVALUES  SPECIFIED  BY  THE  PARAMETER  IORD.  THE 
SUBROUTINE  THEN  SOLVES  THE  LINEAR  SYSTEM 

X  •  E  •  Z11  =  Z21 

FOR  X,  WHERE  Z11  AND  Z21  ARE  THE  UPPER  AND  LOWER  LEFT  N  PY  N 
BLOCKS  OF  Z. 

REF.:  ARNOLD,  W.F.,  "ON  THE  NUMERICAL  SOLUTION  OF 
ALGEBRAIC  MATRIX  RICCATI  EQUATIONS,"  PHD  THESIS,  DSC, 
DECEMBER  1983. 

•••••PARAMETER  DESCRIPTION: 

ON  INPUT: 
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NR  INTEGER 

ROW  DIMENSION  OF  THE  ARRAY  CONTAINING  THE  MATRIX  E 
AS  DECLARED  IN  THE  MAIN  CALLING  PROGRAM  DIMENSION 
STATEMENT; 

NRD  INTEGER 

ROW  DIMENSION  OF  THE  ARRAYS  CONTAINING  THE  MATRICES 
G,  F  AND  Z  AS  DECLARED  IN  THE  MAIN  CALLING  PROGRAM 
DIMENSION  STATEMENT.; 


NN  INTEGER 

ORDER  OF  THE  SQUARE  MATRICES  G  AND  F; 


N  INTEGER 

ORDER  OF  THE  SQUARE  MATRIX  E; 


G  REAL (NRD, NN) 

PENCIL  MATRIX  CORRESPONDING  TO  THE  GENERALIZED  RICCATI 
PROBLEM,  WHOSE  CONTENTS  ARE  ALTERED  BY  THIS  ROUTINE; 


F  REAL(NRD.NN) 

PENCIL  MATRIX  CORRESPONDING  TO  THE  GENERALIZED  RICCATI 
PROBLEM.  WHOSE  CONTENTS  ARE  ALTERED  BY  THIS  ROUTINE; 


E  REAL (NR, N) 

DESCRIPTOR  MATRIX  OR  CO-ORDINATE  BALANCING  MATRIX 
AS  SPECIFIED  BY  THE  PARAMETER  IBAL; 

CPERM  REAL(NN) 

WORKING  VECTOR  OF  SIZE  AT  LEAST  NN; 


CSCALE  REAL(NN) 

WORKING  VECTOR  OF  SIZE  AT  LEAST  NN; 


IND  INTEGER (NN) 

WORKING  VECTOR  OF  SIZE  AT  LEAST  NN; 


IORD  INTEGER 

PARAMETER  SPECIFYING  THE  SPECTRUM  OF  THE  UPPER'  LEFT  N 
BY  N  BLOCK  OF  THE  ORDERED  REAL  SCHUR  FORM  AS  FOLLOWS: 
=  1  GENERALIZED  EIGENVALUES  WHOSE  MAGNITUDE  IS  LESS 

THAN  UNITY 
=  0  ANY  ORDER 

=  -1  GENERALIZED  EIGENVALUES  WHOSE  REAL  PARTS  ARE 
LESS  THAN  ZERO; 

IBAL  INTEGER 

PARAMETER  SPECIFYING  THE  BALANCING  BEING  EMPLOYED  AS 
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FOLLOWS: 

=  0  WARD  BALANCING  AND  E  IS  A  DESCRIPTOR  MATRIX 
=  1  CO-ORDINATE  BALANCING  AND  E  IS  THE  BALANCING 
TRANSFORMATION 

=  2  NO  BALANCING  AND  E  IS  A  DESCRIPTOR  MATRIX; 


TYPE 


LOGICAL 

s  .TRUE.  FOR  CONTINUOUS-TIME  SYSTEM 
=  .FALSE.  FOR  DISCRETE -TIME  SYSTEM; 


EFLAG 


CHARACTER 

FLAG  SET  TO  ’Y’  IF  E  IS  A  DESCRIPTOR  MATRIX  THAT  IS 
OTHER  THAN  THE  IDENTITY  MATRIX. 


ON  OUTPUT: 


ALFR 


ALFI 


BETA 


CONTAINS  THE  SOLUTION  MATRIX  X  COMPUTED  AS  SHOWN 
ABOVE; 


REAL(NRD.NN) 

CONTAINS  THE  MATRIX  PRODUCT 


c 

(  E 

0  ) 

( 

Z11 

Z12) 

c 

( 

) 

•  ( 

) 

c 

(  0 

I  ) 

( 

Z21 

Z22) 

WHERE  Z  IS  THE  ORTHOGONAL  TRANSFORMATION  MATRIX 
DESCRIBED  ABOVE; 


REAL(NN) 

REAL  PARTS  OF  THE  DIAGONAL  ELEMENTS  THAT  WOULD  RESULT 
IF  THE  Q  AND  Z  TRANSFORMATIONS  WERE  APPLIED  TO  THE 
G  MATRIX  SUCH  THAT  IT  WOULD  BE  REDUCED  COMPLETELY  TO 
TRIANGULAR  FORM  AND  THE  DIAGONAL  ELEMENTS  OF  THE 
TRANSFORMED  F  MATRIX  (ALSO  TRIANGULAR)  WOULD  BE  REAL 
AND  POSITIVE; 


REAL(NN) 

IMAGINARY  PARTS  OF  THE  DIAGONAL  ELEMENTS  THAT  WOULD 
RESULT  IF  THE  Q  AND  THE  Z  TRANSFORMATIONS  WERE 
APPLIED  TO  THE  G  MATRIX  SUCH  THAT  IT  WOULD  BE  REDUCED 
COMPLETELY  TO  TRIANGULAR  FORM  AND  THE  DIAGONAL, 
ELEMENTS  OF  THE  TRANSFORMED  F  MATRIX  (ALSO  TRIANGULAR) 
WOULD  BE  REAL  AND  POSITIVE.  NONZERO  VALUES  OCCUR  IN 
PAIRS;  THE  FIRST  MEMBER  IS  POSITIVE  AND  THE  SECOND 
MEMBER  IS  NEGATIVE; 


REAL(NN) 

REAL  NONNEGATIVE  DIAGONAL  ELEMENTS  OF  F  THAT  WOULD 
RESULT  IF  G  WERE  REDUCED  COMPLETELY  TO  TRAINGULAR 
FORM;  THE  GENERALIZED  EIGENVALUES  ARE  THEN  GIVEN  BY 
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THE  RATIOS  (( ALFR  ♦  I«ALFI)/BETA ) ; 


CPERM(I) 

CONDITION  ESTIMATE  OF  E»Z11  WITH  RESPECT  TO  INVERSION 

IND(1)  ERROR  FLAG  AS  FOLLOWS 

=  0  INDICATES  NORMAL  RETURN 

=  NONZERO  IF  MORE  THAT  50  ITERATIONS  WERE  REQUIRED  TO 
DETERMINE  THE  DIAGONAL  BLOCKS  FOR  THE  QUASITRIANGULAR 
FORM; 

IND(2)  ERROR  FUG  AS  FOLLOWS 

=  0  INDICATES  NORMAL  RETURN 
s  1  INDICATES  ATTEMPTED  REORDERING  FAILED. 

••»•• ALGORITHM  NOTES: 

NONE. 


•••••HISTORY: 

THIS  SUBROUTINE  WAS  WRITTEN  BY  W.F.  ARNOLD,  NAVAL  WEAPONS 
CENTER.  CODE  35104,  CHINA  UKE,  CA  93555,  AS  PART  OF  THE 
SOFTWARE  PACKAGE  RICPACK.  SEPTEMBER  1983. 
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SUBROUTINE  NEWT(NR,NRR,NRW,NRX,N,M,E,A,B,CQC,R,S,RI,RSDM,X,W1, 
X  W2 , W3 , ALFR , ALFI , WK1 , IPVT , RTOL , EFLAG , RFLAG , 

X  RDFLG , SFLAG , SE  P , TYPE , MAXIT , INFO , NOUT ) 


•••••PARAMETERS: 

INTEGER  NR,NRR,NRW,NRX,N,M,IPVT(N) , MAXIT, INFO, NOUT 
CHARACTER  EFLAG, RFLAG, RDFLG, SFLAG 

DOUBLE  PRECISION  E(NR.N) ,A(NR,N) ,B(NR,M) ,CQC(NR,N) ,R(NR,M) , 
X  S(NR,M) ,RI(NR,M) ,RSDM(NRR,N) ,X(NRX,N) , 

X  W1(NRW,N),W2(NRW,N),W3(NRW,N),ALFR(N), 

X  ALFI(N),WK1 (N), RTOL, SEP 

LOGICAL  TYPE 

•••••LOCAL  VARIABLES: 

INTEGER  I , IER 1 , IER2 , ITER , J 

DOUBLE  PRECISION  DPN1 ,EPS,E1N,H1N,TOL,T1 .T2.T3.X1N 


•••••FORTRAN  FUNCTIONS: 
DOUBLE  PRECISION  DSQRT 


•••••FUNCTION  SUBPROGRAMS: 

DOUBLE  PRECISION  D1NRM 

•••••SUBROUTINES  CALLED: 

FBGAIN,  LYPCND,  LYPDSD,  MADD,  MMUL,  MQF,  MSCALE,  MSUB,  MULA, 
SEPEST,  TRNATB 


•••••PURPOSE: 

GIVEN  THE  MODEL  MATRICES  FOR  THE  CONTINUOUS-  OR  DISCRETE-TIME 
OPTIMAL  REGULATOR  PROBLEM  THAT  RESULTS  IN  A  GENERALIZED 
ALGEBRAIC  RICCATI  EQUATION  (GARE),  AND  AN  INITIAL  GUESS  FOR 
THE  SOLUTION  TO  THE  GARE,  THIS  SUBROUTINE  APPLIES  A  NEWTON 
TYPE  ITERATIVE  REFINEMENT  PROCEDURE.  TO  GUARANTEE 
CONVERGENCE,  THE  INITIAL  GUESS  MUST  STABILIZE  THE  CLOSED  LOOP 
SYSTEM  MATRIX  E»»-1«(A  -  B»K),  WHERE 


K  =  (R»*-1 )*( BT*X*E  ♦  ST) 


CONTINUOUS 


K  =  ((R  +  BT»X»B)”-1  )»(BT»X»A  +  ST)  DISCRETE 

AND  T  DENOTES  MATRIX  TRANSPOSE. 

REF.:  ARNOLD,  W.F.,  "ON  THE  NUMERICAL  SOLUTION  OF 
ALGEBRAIC  MATRIX  RICCATI  EQUATIONS,"  PHD  THESIS,  USC, 
DECEMBER  1983. 

•••••PARAMETER  DESCRIPTION: 
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ON  INPUT: 


INTEGER 

ROW  DIMENSION  OF  THE  ARRAYS  CONTAINING  THE  MATRICES 
E,  A,  B,  CQC,  R,  S  AND  RI  AS  DECLARED  IN  THE  MAIN 
CALLING  PROGRAM  DIMENSION  STATEMENT; 

INTEGER 

ROW  DIMENSION  OF  THE  ARRAY  CONTAINING  THE  MATRIX 
RSDM  AS  DECLARED  IN'* THE  MAIN  CALLING  PROGRAM 
DIMENSION  STATEMENT; 

INTEGER 

ROW  DIMENSION  OF  THE  ARRAYS  CONTAINING  THE  MATRICES 
W1,  W2  AND  W3  AS  DECLARED  IN  THE  MAIN  CALLING 
PROGRAM  DIMENSION  STATEMENT; 

INTEGER 

ROW  DIMENSION  OF  THE  ARRAY  CONTAINING  THE  MATRIX  X 
AS  DECLARED  IN  THE  MAIN  CALLING  PROGRAM  DIMENSION 
STATEMENT; 

INTEGER 

ORDER  OF  THE  SQUARE  MATRICES  E.  A,  CQC,  RSDM  AND  X 
ROW  DIMENSION  OF  THE  MATRICES  B  AND  S; 

INTEGER 

ORDER  OF  THE  SQUARE  MATRICES  R  AND  RI 
COLUMN  DIMENSION  OF  THE  MATRICES  B  AND  S; 

REAL(NR.N) 

MODEL  DESCRIPTOR  MATRIX; 

REAL (NR, N> 

MODEL  SYSTEM  MATRIX; 

REAL (NR, M) 

MODEL  INPUT  MATRIX; 

REAL (NR, N) 

MATRIX  PRODUCT  CT«Q«C  WHERE  T  DENOTES  MATRIX 
TRANSPOSE ; 


REAL (NR, M) 

CONTROL  WEIGHTING  MATRIX; 

REAL(NR.M) 

STATE  -  INPUT  CROSS-WEIGHTING  MATRIX; 
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REAL (NR, M) 

INVERSE  OF  THE  CONTROL  WEIGHTING  MATRIX; 
REAL(NRX.N) 

INITIAL  GUESS  FOR  RICCATI  SOLUTION  THAT  MUST 
STABILIZE  THE  CLOSED  LOOP  SYSTFM  MATRIX; 


W1,W2,W3 

REAL(NRW.N) 

SCRATCH  ARRAYS  OF  SIZE  AT  LEAST  N  BY  N; 

WK1  REAL(N) 

WORK  VECTOR  OF  LENGTH  AT  LEAST  N; 

I PVT  INTEGER (N) 

WORK  VECTOR  OF  LENGTH  AT  LEAST  N; 

RTOL  REAL 

TOLERANCE  ON  THE  CONDITION  ESTIMATE  OF  R+BT»X*B 
WITH  RESPECT  TO  INVERSION  (DISCRETE  PROBLEM). 

ERROR  RETURN  IF  THIS  TOLERANCE  IS  EXCEEDED; 

EFLAG  CHARACTER 

FLAG  SET  TO  ’Y'  IF  E  IS  OTHER  THAN  THE  IDENTITY 
MATRIX; 

RFLAG  CHARACTER 

FLAG  SET  TO  'Y'  IF  R  IS  OTHER  THAN  THE  IDENTITY 
MATRIX; 

RDFLG  CHARACTER 

FLAG  SET  TO  'Y'  IF  R  IS  A  DIAGONAL  MATRIX; 

SFLAG  CHARACTER 

FLAG  SET  TO  'Y'  IF  S  IS  OTHER  THAN  THE  ZERO  MATRIX; 
TYPE  LOGICAL 

=  .TRUE.  FOR  CONTINUOUS-TIME  SYSTEM 
=  .FALSE.  FOR  DISCRETE-TIME  SYSTEM; 

MAXIT  INTEGER 

MAXIMUM  NUMBER  OF  NEWTON  ITERATIONS  PERMITTED; 

NOUT  INTEGER 

UNIT  NUMBER  OF  OUTPUT  DEVICE  FOR  ERROR  WARNING 
MESSAGES. 

ON  OUTPUT: 
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RSDM  REAL (NR R, N) 

DIFFERENCE  MATRIX  BETWEEN  SOLUTIONS  AT  LAST  TWO 
ITERATIONS ; 

X  REFINED  RICCATI  SOLUTION  MATRIX; 

ALFR  REAL(N) 

REAL  PARTS  OF  THE  CLOSED  LOOP  SYSTEM  EIGENVALUES 

UNDER  THE  OPTIMAL  FEEDBACK  FOR  THE  RICCATI  SOLUTION 
AT  THE  LAST  ITERATION; 

ALFI  REAL(N) 

IMAGINARY  PARTS  OF  THE  CLOSED  LOOP  SYSTEM 
EIGENVALUES  UNDER  THE  OPTIMAL  FEEDBACK  FOR  THE 
RICCATI  SOLUTION  AT  THE  LAST  ITERATION; 

WK1 (1 )  1-NORM  OF  THE  RSDM  MATRIX; 

IPVT(1)  CONVERGENCE  CRITERIA  INDICATOR; 

SEP  REAL 

ESTIMATE  OF  THE  SEPARATION  OF  THE  CLOSED  LOOP 
SPECTRUM  AT  THE  LAST  ITERATION  (CONTINUOUS  PROBLEM); 

MAXIT  NUMBER  OF  ITERATIONS  PERFORMED; 

INFO  INTEGER 

ERROR  FLAG  WITH  THE  FOLLOWING  MEANINGS 
=  -1  NO  CONVERGENCE 
=  -2  INITIAL  GUESS  NOT  STABILIZING 
=  -3  E  MATRIX  NOT  IDENTITY 

s  -4  INDICATES  A  FAILURE  OF  THE  QR  ALGORITHM  TO 
DETERMINE  THE  EIGENVALUES  IN  SOLVING  THE 
LYAPUNOV  EQUATION 

=  -5  CONDITION  ESTIMATE  OF  R+BT»X»B  WITH  RESPECT  TO 
INVERSION  EXCEEDS  THE  TOLERANCE  VALUE  RTOL. 

•••••ALGORITHM  NOTES: 

THE  ALGORITHM  CURRENTLY  EMPLOYED  IS  BASED  ON  THE  BARTELS  ►- 
STEWART  ALGORITHM  FOR  LYAPUNOV  EQUATIONS  AND  IS  VALID  FOR  THE 
CASE  E  =  IDENTITY  ONLY  AT  THIS  TIME. 

•••••HISTORY: 

THIS  SUBROUTINE  WAS  WRITTEN  BY  W.F.  ARNOLD,  NAVAL  WEAPONS 
CENTER,  CODE  35104,  CHINA  LAKE,  CA  93555,  AS  PART  OF  THE 
SOFTWARE  PACKAGE  RICPACK,  SEPTEMBER  1983. 
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SUBROUTINE  ORDER  (A,B,Z,NMAX,N .EPS, IFAIL, TYPE , IFIRST, IND) 


•••••PARAMETERS: 

INTEGER  NMAX.N, IFAIL, IFIRST, IND (N) 

DOUBLE  PRECISION  A(NMAX,N) ,B(NMAX,N) ,Z(NMAX,N) .EPS 
LOGICAL  TYPE 


•••••LOCAL  VARIABLES: 

INTEGER  I, II, III , IS, ISTEP ,K,L,LS , LS 1 , LS2,L1 ,L2 ,NUM 


•••••FORTRAN  FUNCTIONS: 
DOUBLE  PRECISION  DABS 


•••••SUBROUTINES  CALLED: 
EXCHQZ 


»••*• PUR POSE: 

GIVEN  THE  UPPER  TRIANGULAR  MATRIX  B  AND  UPPER  HESSENBERG 
MATRIX  A  WITH  1X1  OR  2X2  DIAGONAL  BLOCKS,  THIS  SUBROUTINE 
REORDERS  THE  DIAGONAL  BLOCKS  ALONG  WITH  THE  GENERALIZED 
EIGENVALUES  CORRESPONDING  TO  THE  REGULAR  MATRIX  PENCIL 
A  -  LAMBDA »B  BY  CONSTRUCTING  ROW  AND  COLUMN  EQUIVALENCE 
TRANSFORMATIONS  QT  AND  ZT.  THE  COLUMN  TRANSFORMATIONS  ARE 
THEN  APPLIED  TO  THE  MATRIX  Z. 

REF.:  VAN  DOOREN,  P.,  A  GENERALIZED  EIGENVALUE  APPROACH  FOR 
SOLVING  RICCATI  EQUATIONS,  SIAM  J.  SCI.  STAT.  COMPUT., 

VOL.  2,  NO.  2,  JUNE  1981,  121-135. 


•PARAMETER  DESCRIPTION: 


ON  INPUT: 


INTEGER 

ROW  DIMENSION  OF  THE  ARRAYS  CONTAINING  A,B,Z  AS 
DECLARED  IN  THE  MAIN  CALLING  PROGRAM  DIMENSION 
STATEMENT; 

INTEGER 

ORDER  OF  THE  MATRICES  A,B,Z; 

REAL(NMAX.N) 

UPPER  HESSENBERG  MATRIX  WITH  1X1  OR  2X2  DIAGONAL 
BLOCKS.  ELEMENTS  OUTSIDE  THE  UPPER  HESSENBERG 
STRUCTURE  ARE  ARBITRARY; 

REAL (NMAX.N) 

UPPER  TRIANGULAR  MATRIX.  ELEMENTS  OUTSIDE  THE 
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UPPER  TRIANGULAR  STRUCTURE  ARE  ARBITRARY; 

EPS  REAL 

REQUIRED  ABSOLUTE  ACCURACY  OF  THE  RESULTS. 

NORMALLY  EQUAL  TO  THE  MACHINE  PRECISION; 

TYPE  LOGICAL 

CONTROL  PARAMETER  THAT  SPECIFIES  THE  REGIONS  OF  THE 
THE  COMPLEX  PLANE  THAT  THE  GENERALIZED  EIGENVALUES 
ARE  ORDERED  BY.  TO  CONTROL  THE  REGION  THAT  APPEARS 
FIRST,  SEE  IFIRST  BELOW 

=  .TRUE.  GENERALIZED  EIGENVALUES  ARE  ORDERED  BY 
REGION  INSIDE  THE  COMPLEX  LEFT  HALF  PLANE  OR 
OUTSIDE  THIS  REGION 

=  .FALSE.  GENERALIZED  EIGENVALUES  ORDERED  BY  REGION 
INSIDE  THE  UNIT  CIRCLE  OR  OUTSIDE  THIS  REGION; 

IFIRST  INTEGER 

CONTROL  PARAMETER  THAT  SPECIFIES  WHICH  OF  THE 
REGIONS  SPECIFIED  BY  TYPE (SEE  ABOVE)  APPEARS  FIRST 
(I.E.  IN  THE  UPPER  LEFT  NXN  BLOCK) 

=  -1  INSIDE  REGION  APPEARS  FIRST 
=  +1  OUTSIDE  REGION  APPEARS  FIRST 
IFIRST  IS  ALTERED  BY  THIS  SUBROUTINE; 

IND  INTEGER <N) 

WORKING  ARRAY  THAT  IS  ALTERED  BY  THIS  SUBROUTINE. 


ON  OUTPUT: 

A.B  UPPER  HESSENBERG  MATRIX,  UPPER  TRIANGULAR  MATRIX 

REORDERED  AS  SPECIFIED  BY  TYPE  AND  IFIRST (ABOVE ) ; 

Z  REAL(NMAX.N) 

THIS  ARRAY  IS  OVERWRITTEN  BY  THE  PRODUCT  OF  THE 
CONTENTS  OF  THE  ARRAY  Z(UPON  ENTRY  INTO  THIS 
SUBROUTINE),  AND  THE  COLUMN  TRANSFORMATIONS  ZT 
(CALCULATED  BY  THIS  SUBROUTINE); 

IFAIL  INTEGER 

ERROR  FLAG 

=  1  INDICATES  ATTEMPTED  REORDERING  FAILED 
=  0  NORMAL  RETURN. 

*****  ALGORITHM  NOTES: 

NONE. 


•••••HISTORY: 

ORIGINAL  VERSION  THAT  SORTED  BY  UNIT  CIRCLE  REGION  OF  COMPLEX 
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PLANE  WRITTEN  BY  P.  VAN  DOOREN ("A  GENERALIZED  EIGENVALUE 
APPROACH  FOR  SOLVING  RICCATI  EQUATIONS",  INTERNAL  REPORT 
NA -80-02 ,  DEPT.  OF  COMPUTER  SCIENCE,  STANFORD  UNIVERSITY, 
1980).  THIS  VERSION  MODIFIED  BY  W.  F.  ARNOLD (DEPT.  OF 
ELECTRICAL  ENGINEERING  -  SYSTEMS,  UNIV.  OF  SOUTHERN  CALIF 
LOS  ANGELES,  CA  90089)  TO  INCLUDE  THE  SORTING  CONTROL 
PARAMETER  "TYPE",  SEPT  1982. 


SUBROUTINE  LYPCND ( NF , NH , N , F , H . Z , MR , WI , WK , IER 1 , IER2 ) 


•••••PARAMETERS: 

INTEGER  NF ,  NH ,  N , IER 1 , IER  2 

DOUBLE  PRECISION  F(NF,N) ,H(NH,N) ,Z(NF,N) ,WR(N) ,WI(N) ,WK(N) 

•••••LOCAL  VARIABLES: 

INTEGER  LOW.IGH 

•••••subrouhnes  called: 

ORTHES , ORTRAN , HQRORT , MQFWO (MULWOA ) .SYMSLV (LINEQ , DGECOM , DGESLM) 
TRNATA 


•••••PURPOSE: 

THIS  SUBROUTINE  SOLVES  THE  CONTINUOUS  TIME  LYAPUNOV  EQUATION 
T 

F  »X  +  X»F  ♦  H  s  0. 

BY  THE  BARTELS-STEWART  ALGORITHM  (SEE  REF.(1)). 

•••••PARAMETER  DESCRIPTION: 

ON  INPUT: 

NF.NH  ROW  DIMENSIONS  OF  THE  ARRAYS  CONTAINING  F 

(AND  Z)  ,H,  RESPECTIVELY,  AS  DECLARED  IN 
MAIN  CALLING  PROGRAM  DIMENSION  STATEMENT; 

N  ORDER  OF  THE  MATRICES  F  AND  H; 

F  N  X  N  (REAL)  MATRIX; 

H  N  X  N  SYMMETRIC  MATRIX; 

IER1  INTEGER  VARIABLE;  NORMALLY  SET  IER1  TO  0; 

IF  IER 1  IS  SET  TO  A  NON-ZERO  INTEGER,  THE 
REDUCTION  OF  F  TO  REAL  SCHUR  FORM  IS ■ 
SKIPPED  AND  THE  ARRAYS  F  AND  Z  ARE  ASSUMED 
TO  CONTAIN  THE  REAL  SCHUR  FORM  AND 
ACCOMPANYING  ORTHOGONAL  MATRIX  THUS 
PERMITTING  MORE  EFFICIENT  SOLUTION  OF 
SEVERAL  EQUATIONS  WITH  DIFFERENT  CONSTANT 
TERMS  H; 


IER2 


INTEGER  VARIABLE;  NORMALLY  SET  IER2  TO  0 
IF  ONLY  A  REAL  SCHUR  FORM  OF  F  AND 
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ASSOCIATED  ORTHOGONAL  SIMILARITY  MATRIX  Z 
ARE  DESIRED  SET  IER2  TO  A  NON-ZERO 
INTEGER. 

ON  OUTPUT: 

H  N  X  N  ARRAY  CONTAINING  THE  (SYMMETRIC) 

SOLUTION  X  OF  THE  LYAPUNOV  EQUATION; 

F  N  X  N  ARRAY  CONTAINING  IN  ITS  UPPER 

TRIANGLE  AND  FIRST  SUBDIAGONAL  A  REAL 
SCHUR  FORM  OF  F; 

Z  N  X  N  ARRAY  CONTAINING,  ON  OUTPUT,  THE 

ORTHOGONAL  MATRIX  THAT  REDUCES  F  TO  REAL 
SCHUR  FORM; 

WR  REAL  SCRATCH  VECTOR  OF  LENGTH  N;  ON  OUTPUT 

(WR(I),I*1,N)  CONTAINS  THE  REAL  PARTS  OF 
THE  EIGENVALUES  OF  F  AND  THUS  CAN  BE  USED 
TO  TEST  THE  STABILITY  OF  F; 

WI  REAL  SCRATCH  VECTOR  OF  LENGTH  N;  ON  OUTPUT 

CONTAINS  THE  IMAGINARY  PARTS  OF  THE 
EIGENVALUES  OF  F; 

WK  REAL  SCRATCH  VECTOR  OF  LENGTH  N; 

IER1  =0  FOR  NORMAL  RETURN  (IF  =0  ON  INPUT), 

=J  IF  THE  J-TH  EIGENVALUE  HAS  NOT  BEEN 
DETERMINED  IN  THE  QR  ALGORITHM  (IF  =0  ON 
INPUT). 

•••••ALGORITHM  NOTES: 

IT  IS  ASSUMED  THAT  F  HAS  NO  EIGENVALUES  WHICH  SUM  TO  ZERO  (THIS 

CAN  BE  CHECKED  FROM  THE  ARRAY  WR).  THIS  IS  SUFFICIENT  TO 

GUARANTEE  A  UNIQUE  SOLUTION. 

IF,  MOREOVER,  F  IS  STABLE  THEN.  X  IS  NONNEGATIVE  DEFINITE. 

*«»t*R£f-gR£NC£S  : 

(1)  BARTELS,  R.H.,  AND  G.W.  STEWART,  SOLUTION 
OF  THE  MATRIX  EQUATION  AX  +  XB  =  C, 

ALGORITHM  432,  COMM.  ACM,  15(1 972) ,820-826 . 


•••••HISTORY: 

WRITTEN  BY  ALAN  J.  LAUB  (DEP'T.  OF  EE-SYSTEMS,  U.  OF  SOUTHERN 
CALIF.,  LOS  ANGELES,  CA  90089,  PH.:  (213)  743-5535)  SEP.  1977 
MOST  RECENT  VERSION:  JUNE  28,  1982. 
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SUBROUTINE  LYPDSD (NF , NH , N , F , H , Z , WR , WI , WK , U , IDIM , IER 1 , IER2 ) 
•••••PARAMETERS: 

INTEGER  NF,NH,N,IER1,IER2,IDIM(N) 

DOUBLE  PRECISION  F(NF,N) ,U(NF,N) ,H(NH,N) ,Z(NF,N) ,WR(N) ,WI(N) , 

X  WK(N) 

•••••LOCAL  VARIABLES: 

INTEGER  LOW , IGH , KIN , KOUT 

•••••SUBROUTINES  CALLED: 

ORTHES , ORTRAN , HQRORT , MQFWO ( MULWOA ) , DSTSLV (LINEQ , DDCOM P , DSOLVE ) 
TRNATA 


•••••PURPOSE: 
THIS  SUBROUTINE 

SOLVES  THE  DISCRETE  TIME  LYAPUNOV  EQUATION 

T 

F  »X»F-  X 

=  H. 

BY  A  MODIFICATION  OF  THE  BARTELS-STEWART  ALGORITHM  (SEE  REFS. 
(1)  AND  (2)). 

•••••PARAMETER 
ON  INPUT: 

DESCRIPTION: 

(AND  Z, U) .AND  H  .RESPECTIVELY,  AS  DECLARED 
IN  THE  CALLING  PROGRAM  DIMENSION  STATEMENT 

N 

ORDER  OF  THE  MATRICES  F  AND  C; 

F 

N  X  N  (REAL)  MATRIX; 

H 

N  X  N  SYMMETRIC  MATRIX; 

IER1 

INTEGER  VARIABLE;  NORMALLY  SET  IER1  TO  0; 
IF  IER 1  IS' SET  TO  A  NON-ZERO  INTEGER,  THE 
REDUCTION  OF  F  TO  REAL  SCHUR  FORM  IS, 
SKIPPED  AND  THE  ARRAYS  F  AND  Z  ARE  ASSUMED 
TO  CONTAIN  THE  REAL  SCHUR  FORM  AND 
ACCOMPANYING  ORTHOGONAL  MATRIX  THUS 
PERMITTING  MORE  EFFICIENT  SOLUTION  OF 
SEVERAL  EQUATIONS  WITH  DIFFERENT  CONSTANT 
TERMS  H; 

IER2 

INTEGER  VARIABLE;  NORMALLY  SET  IER2  TO  0; 
IF  ONLY  A  REAL  SCHUR  FORM  OF  F  AND 

m 


OH  OUTPUT: 
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ASSOCIATED  ORTHOGONAL  SIMILARITY  MATRIX  Z 
ARE  DESIRED  SET  IER2  TO  A  NON-ZERO 
INTEGER. 


N  X  N  ARRAY  CONTAINING  THE  (SYMMETRIC) 
SOLUTION  X  OF  THE  LYAPUNOV  EQUATION; 


N  X  N  ARRAY  CONTAINING  IN  ITS  UPPER 
TRIANGLE  AND  FIRST  SUBDIAGONAL  A  REAL 
SCHUR  FORM  OF  F; 


N  X  N  ARRAY  CONTAINING,  ON  OUTPUT,  THE 
ORTHOGONAL  MATRIX  THAT  REDUCES  F  TO  REAL 
SCHUR  FORM; 


REAL  SCRATCH  VECTOR  OF  LENGTH  N;  ON  OUTPUT 
(WR(I) ,1=1 ,N)  CONTAINS  THE  REAL  PARTS  OF 
THE  EIGENVALUES  OF  F  AND  THUS  CAN  BE  USED 
TO  TEST  THE  STABILITY  OF  F; 


REAL  SCRATCH  VECTOR  OF  LENGTH  N;  ON  OUTPUT 
CONTAINS  THE  IMAGINARY  PART  OF  THE 
EIGENVALUES  OF  F; 


REAL  SCRATCH  VECTOR  OF  LENGTH  N; 


N  X  N  REAL  SCRATCH  ARRAY; 


INTEGER  SCRATCH  VECTOR  OF  LENGTH  N; 


=0  FOR  NORMAL  RETURN  (IF  =0  ON  INPUT), 

=J  IF  THE  J-TH  EIGENVALUE  HAS  NOT  BEEN 
DETERMINED  IN  THE  QR  ALGORITHM  (IF  =0  ON 
INPUT). 


•••••ALGORITHM  NOTES: 

IT  IS  ASSUMED  THAT  F  HAS  NO  EIGENVALUES  WITH  PRODUCT  EQUAL  TO 
ZERO  (THIS  CAN  BE  CHECKED  FROM  THE  ARRAY  WR).  THIS  IS 
SUFFICIENT  TO  GUARANTEE  A  UNIQUE  SOLUTION. 

IF,  MOREOVER,  F  IS  STABLE  THEN  X  IS  NONNEGATIVE  DEFINITE. 


'••REFERENCES: 

(1)  BARTELS,  R.H.,  AND  G'.W.  STEWART,  SOLUTION 
OF  THE  MATRIX  EQUATION  AX  ♦  XB  =  C, 
ALGORITHM  432,  COMM.  ACM,  15(1972) ,820-826, 


(2)  BARRAUD.A. Y. ,  A  NUMERICAL  ALGORITHM  TO  SOLVE  A  XA-X=Q, 


oooonooo 
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IEEE  TRANSACTIONS  ON  AUTOMATIC  CONTROL,  VOL.  AC-22, 
NO. 5,  OCTOBER  1977,883-885. 


•••••HISTORY: 

WRITTEN  BY  J.A.K.  CARRIG  (ELEC.  SYS.  LAB.,  M.I.T.,  RM.  35-427, 
CAMBRIDGE,  MA  02139,  PH.:  (617)  653-7263,  SEPTEMBER  1978. 

MOST  RECENT  VERSION:  SEPT.  20,  1978. 


on  o  o  oo  n  o 


NWC  TP  6521 


SUBROUTINE  SEPEST ( NR , N , T , Q , SE P , INFO ) 

•••••PARAMETERS: 

INTEGER  NR, N, INFO 

DOUBLE  PRECISION  T(NR,N) ,Q(NR,N) ,SEP 
•••••LOCAL  VARIABLES: 

INTEGER  IPVT(4) ,I,II,IM1 ,J, JP1 ,K1 .K2.K1M1.L1 .L2.L1M1 ,L1MK1 , 
X  M.ND.NM1 

DOUBLE  PRECISION  A(4,4)  ,VEC(4.>,Z(4)  ,A1N,RCOND,S,TEMP,T1N 

•••••FORTRAN  FUNCTIONS: 

DOUBLE  PRECISION  DABS , DMAX1 , DMIN 1 , DSIGN 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


•••••FUNCTION  SUBPROGRAMS: 
DOUBLE  PRECISION  D1NRM 

•••••SUBROUTINES  CALLED: 
DGECOM,  DGESLM,  MSCALE,  SYMSLV 


•••••PURPOSE: 

GIVEN  A  QUASITRIANGULAR  MATRIX  T  AND  A  SYMMETRIC  MATRIX  Q, 

THIS  SUBROUTINE  COMPUTES 

SEP  =  MIN (  1-NORM(  T -TRANSPOSED  ♦  Q*T)/  l-NORM(Q)) 

REF.:  BARTELS,  R.H.  AND  G.W.  STEWART,  "SOLUTION  OF  THE  MATRIX 
EQUATION  A»X  ♦  X»B  =  C,"  COMM.  OF  THE  ACM,  VOL.  15, 

PP.  820-826,  1972. 

CLINE,  A.K.,  MOLER,  C.B.,  STEWART,  G.W.  AND  J.H.  WILKINSON, 

"AN  ESTIMATE  OF  THE  CONDITION  NUMBER  OF  A  MATRIX,"  SIAM  J.  OF 
NUMERICAL  ANALYSIS,  VOL.  16,  PP.  368-375,  1979. 

•••••PARAMETER  DESCRIPTION: 

ON  INPUT: 

NR  INTEGER 

ROW  DIMENSION  OF  THE  ARRAYS  CONTAINING  THE  MATRICES 
T  AND  Q  AS  DECLARED  IN  THE  MAIN  CALLING  PROGRAM 
DIMENSION  STATEMENT; 

N  INTEGER 

ORDER  OF  THE  SQUARE  MATRICES  Q  AND  T; 

T  REAL (NR, N) 

QUASITRIANGULAR  INPUT  MATRIX. 
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C 

C  ON  OUTPUT: 

C 

C  SEP  REAL 

C  AN  ESTIMATE  OF  THE  QUANTITY  SPECIFIED  ABOVE; 

C 

C  Q  REAL (NR, N) 

C  THE  MATRIX  THAT  MINIMIZES  THE  QUANTITY  SEP; 

C 

C  INFO  INTEGER 

C  ERROR  FLAG  WITH  THE  FOLLOWING  MEANINGS 

C  =0  INDICATES  NORMAL  RETURN 

C  s  (N+1 )*L+K  INDICATES  THE  L-TH  AND  K-TH  EIGENVALUES 

C  OF  T  FORM  A  +/-  PAIR,  SO  SEP  IS  EQUAL  TO  ZERO. 

C 

C  •••••ALGORITHM  NOTES: 

C  T 

C  LET  PHI (Y)  =  T»Y  +  Y*T  .  Q  IS  OBTAINED  BY  INVERSE  ITERATION  ON 
C  T 

C  PHI*PHI  .  THE  STARTING  VALUE  OF  Q  IS  CHOSEN  AS  FOLLOWS: 

C  PARTITION  ALL  MATRICES  CONFORMALLY  WITH  T.  Q  IS  CHOSEN  TO 

C  SATISFY 

C  T  T 

C  T»Y  ♦  Y*T  =  B,  T  *Q  ♦  Q»T  .=  Y/1-N0RMCY) . 

C 

c  •••••HISTORY: 

C  THIS  SUBROUTINE  IS  A  MODIFIED  VERSION  OF  THE  SUBROUTINE  OF  THE 

C  SAME  NAME  WRITTEN  BY  RALPH  BEYERS,  2/82  REF.:  BEYERS.  R., 

C  "HAMILTONIAN  AND  SYMPLECTIC  ALGORITHMS  FOR  THE  ALGEBRAIC 

C  RICCATI  EQUATION,"  PHD  THESIS,  CORNELL  UNIVERSITY,  PP. 28 9-2 95, 

C  JANUARY  1983.  THE  MODIFICATIONS  WERE  MADE  BY  W.F.  ARNOLD, 

C  NAVAL  WEAPONS  CENTER,  CODE  35104,  CHINA  LAKE,  CA  93555,  AS 
C  PART  OF  THE  SOFTWARE  PACKAGE  RICPACK,  SEPTEMBER  1983. 

C 
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SUBROUTINE  FBGAIN  (NR,NRX,NRW,N,M,A,B,E,R,RI,S,X,FB,W,WK,IPVT, 

EFLAG , RDFLG , RFLAG , SFLAG , TYPE ) 


•••••PARAMETERS: 

INTEGER  NR  t  NRX , NRW , N , M , IPVT ( N ) 

CHARACTER  EFLAG, RDFLG, RFLAG, SFLAG 

DOUBLE  PRECISION  A(NR,N) ,B(NR,M) ,E(NR,N) ,R(NR,M) ,RI(NR,M) , 
S(NR,M) ,X(NRX,N) ,FB(NRW,N) ,W(NRW,N) ,WK(N) 

LOGICAL  TYPE 

•••••LOCAL  VARIABLES: 

INTEGER  I.J 

DOUBLE  PRECISION  COND 

•••••FORTRAN  FUNCTIONS: 

NONE. 


•••••SUBROUTINES  CALLED: 

MADD,  MLINEQ ,  MMUL,  MULA,  MULB,  TRNATA,  TRNATB 


•••••PURPOSE: 

GIVEN  THE  RICCATI  SOLUTION  AND  THE  MODEL  MATRICES  OF  THE 
OPTIMAL  CONTROL  PROBLEM,  THIS  SUBROUTINE  CALCULATES  THE 
OPTIMAL  FEEDBACK  GAIN  MATRIX  FOR  THE  GENERALIZED  CONTINUOUS- 
OR  DISCRETE -TIME  OPTIMAL  CONTROL  PROBLEM. 

CONTINUOUS:  FB  =  RI»(BT»X»E  ♦  ST) 

DISCRETE:  FB  =  ((R  ♦  BT»X»B)»»-1  )»(BT»X»A  «■  ST) 

WHERE  T  DENOTES  THE  MATRIX  TRANSPOSE. 

REF.:  ARNOLD,  W.F.,  "ON  THE  NUMERICAL  SOLUTION  OF 
ALGEBRAIC  MATRIX  RICCATI  EQUATIONS,"  PHD  THESIS,  USC, 

DECEMBER  1983. 

•••••PARAMETER  DESCRIPTION: 

ON  INPUT: 

NR  INTEGER 

ROW  DIMENSION  OF  THE  ARRAYS  CONTAINING 
A,  B,  E,  R,  RI  AND  S  AS  DECLARED  IN  THE  MAIN 
CALLING  PROGRAM  DIMENSION  STATEMENT; 

NRX  INTEGER 

ROW  DIMENSION  OF  THE  ARRAY  CONTAINING  X  AS  DECLARED 
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IN  THE  MAIN  CALLING  PROGRAM  DIMENSION  STATEMENT: 
INTEGER 

ROW  DIMENSION  OF  THE  ARRAYS  CONTAINING  FB  AND  W  AS 
DECLARED  IN  THE  MAIN  PROGRAM  DIMENSION  STATEMENT; 

INTEGER 

ORDER  OF  THE  SQUARE  MATRICES  A,  E,  AND  X 
ROW  DIMENSION  OF  THE  MATRICES  B,  S,  AND  FB; 


INTEGER 

ORDER  OF  THE  SQUARE  MATRICES  R  AND  RI 
COLUMN  DIMENSION  OF  THE  MATRICES  B  AND  S; 

REAL (NR. N) 

MODEL  SYSTEM  MATRIX; 

REAL (NR, M) 

MODEL  INPUT  MATRIX; 


REAL (NR, N) 

MODEL  DESCRIPTOR  MATRIX; 
REAL (NR. M) 

INPUT  WEIGHTING  MATRIX; 


REAL(NR.M) 

INVERSE  OF  THE  INPUT  WEIGHTING  MATRIX; 

REAL (NR, M) 

STATE  -  INPUT  CROSS-WEIGHTING  MATRIX; 
REAL(NRX.N) 

ALGEBRAIC  RICCATI  EQUATION  SOLUTION  MATRIX; 
REAL(NRW.N) 

SCRATCH  ARRAY  OF  SIZE  AT  LEAST  N  BY  N; 
REAL(N) 

WORKING  VECTOR  OF  LENGTH  AT  LEAST  N; 


IPVT  INTEGER (M) 

WORKING  VECTOR  OF  LENGTH  AT  LEAST  M; 

EFLAG  CHARACTER 

FUG  SET  TO  ’Y’  IF  E  IS  OTHER  THAN  THE  IDENTITY 
MATRIX; 

RDF LG  CHARACTER 


HR 
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C  FUG  SET  TO  'Y'  IF  R  IS  A  DIAGONAL  MATRIX; 

C 

C  RFUG  CHARACTER 

C  FUG  SET  TO  'Y'  IF  R  IS  OTHER  THAN  THE  IDENTITY 

C  MATRIX; 

C 

C  SFUG  CHARACTER 

C  FUG  SET  TO  *Y'  IF  S  IS  OTHER  THAN  THE  ZERO  MATRIX; 

C 

C  TYPE  LOGICAL 

C  =  .TRUE.  FOR  CONTINUOUS-TIME  SYSTEM 

C  =  .FALSE.  FOR  DISCRETE -TIME  SYSTEM. 

C 

C  ON  OUTPUT: 

C 

C  FB  REAL(NRW.N) 

C  OPTIMAL  FEEDBACK  GAIN  MATRIX  AS  DESCRIBED  ABOVE; 

C 

C  WK(1)  ESTIMATED  CONDITION  NUMBER  OF  R+BT*X*B  WITH  RESPECT 

C  TO  INVERSION  (DISCRETE  PROBLEM). 

C 

C  ***** ALGORITHM  NOTES: 

C  NONE. 

C 

c  •••••HISTORY: 

C  THIS  SUBROUTINE  WAS  WRITTEN  BY  W.F.  ARNOLD,  NAVAL  WEAPONS 

C  CENTER,  CODE  35104,  CHINA  UKE,  CA  93555,  AS  PART  OF  THE 

C  SOFTWARE  PACKAGE  RICPACK,  SEPTEMBER  1983. 

C 
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SUBROUTINE  CMPRS(NR,NRDtNRT,N,NN,NNPM,MtE,A,B,CQC,R,S,G,F,U, 
X  MK, WK1 ,WK2 ,WK3 , EFLAG , SFLAG , IBAL , INFO ) 


•••••PARAMETERS: 

INTEGER  NR , NRD , NRT , N , NN , NNPM , M , IBAL , INFO 
CHARACTER  EFLAG, SFLAG 

DOUBLE  PRECISION  £(NR,N) ,A(NR ,N) ,B(NR,M) ,CQC(NR,N) , 

X  R(NR,M) ,S(NR,M) ,G(NRD,NN) ,F(NRD,NN) ,U (NRT, NNPM) , 

X  WK(NRT,M) ,WK1 (M) ,WK2 (M) ,WK3 (NNPM) 

•••••LOCAL  VARIABLES: 

INTEGER  I , J , JOB , K , NN  PI . NPI , NP J , NPK 
C 

C  •••••FORTRAN  FUNCTIONS: 

C  NONE. 

C 

C  •••••SUBROUTINES  CALLED: 

C  DSVDC 

C 

C  - 

c 

C  •••••PURPOSE: 

C  THIS  SUBROUTINE  EMPLOYS  THE  SINGULAR  VALUE  DECOMPOSITION  TO 

C  DETERMINE  AN  ORTHOGONAL  MATRIX  U,  (2N+M)  BY  (2N+M) ,  SUCH  THAT 

C 
C 

C  (  (  )  (  B  )  (  0  ) 

C  (  U11  (  U12  )  (  )  (  ) 

C  (  (  )  (-S  )  =  (  0  ) 

C  ( - ( - )  ( - )  ( - ) 

C  (  U21  (  U22  )  (  R  )  (  RB) 

C 

C  AND  THEN  FORMS  THE  MATRIX  PENCIL 

C 

C  (  E  0  )  (  (A  0  )  ) 

C  LAMBDA *U11*(  )  -  (  U11»(  )  +  U12»(ST  BT)  ) 

C  (  0  ET  )  (  (-CQC  -AT  )  ) 

C 

C  =:  LAMBDA  •  F  -  G 

C 

C  WHERE  T  DENOTES  MATRIX  TRANSPOSE. 

C  THIS  PENCIL  CAN  THEN  BE  USED  FOR  SOLVING  THE  CONTINUOUS-TIME 
C  GARE. 

C 

C  REF.:  ARNOLD,  W.F.,  "ON  THE  NUMERICAL  SOLUTION  OF 
C  ALGEBRAIC  MATRIX  RICCATI  EQUATIONS,"  PHD  THESIS,  USC, 

C  DECEMBER  1983. 

C 

C  •••••PARAMETER  DESCRIPTION: 


n  n 
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ON  INPUT: 


INTEGER 

ROW  DIMENSION  OF  THE  ARRAYS  CONTAINING  THE  MATRICES 
E,  A,  B,  CQC,  R  AND  S  AS  DECLARED  IN  THE  MAIN 
CALLING  PROGRAM  DIMENSION  STATEMENT; 

INTEGER 

ROW  DIMENSION  OF  THE  ARRAYS  CONTAINING  THE  MATRICES 
G  AND  F  AS  DECLARED  IN  THE  MAIN  CALLING  PROGRAM 
DIMENSION  STATEMENT; 

INTEGER 

ROW  DIMENSION  OF  THE  ARRAYS  CONTAINING  THE  MATRICES 
U  AND  WK  AS  DECLARED  IN  THE  MAIN  CALLING  PROGRAM 
DIMENSION  STATEMENT; 

INTEGER 

ORDER  OF  THE  SQUARE  MATRICES  E,  A  AND  CQC 
ROW  DIMENSION  OF  THE  MATRICES  B  AND  S; 

INTEGER 

ORDER  OF  THE  SQUARE  MATRICES  G  AND  F; 

INTEGER 
=  NN  ♦  M; 

INTEGER 

ORDER  OF  THE  SQUARE  MATRIX  R 

COLUMN  DIMENSION  OF  THE  MATRICES  B  AND  S; 

REAL (NR, N) 

MODEL  DESCRIPTOR  MATRIX; 

REAL (NR, N) 

MODEL  SYSTEM  MATRIX; 

REAL (NR, M) 

MODEL  INPUT  MATRIX; 

REAL (NR, N) 

MATRIX  PRODUCT  CT*Q*C  WHERE  T  DENOTES  MATRIX 
TRANSPOSE ; 

REAL(NR.M) 

CONTROL  WEIGHTING  MATRIX; 

REAL(NR,M) 
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STATE  -  INPUT  CROSS-WEIGHTING  MATRIX; 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


WK  REAL(NRT.M) 

SCRATCH  ARRAY  OF  SIZE  AT  LEAST  (NN+M)  BY  M; 

WK1.WK2  REAL(M) 

WORKING  VECTORS  OF  SIZE  AT  LEAST  M; 

WK3  REAL(NNPM) 

WORKING  VECTOR  OF  SIZE  AT  LEAST  NN+M; 

EFLAG  CHARACTER 

FLAG  SET  TO  »Y'  IF  E  IS  OTHER  THAN  THE  IDENTITY 
MATRIX; 

SFLAG  CHARACTER 

FLAG  SET  TO  »Y'  IF  S  IS  OTHER  THAN  THE  ZERO  MATRIX; 
IBAL  INTEGER 

PARAMETER  SET  TO  1  IF  CO-ORDINATE  BALANCING  IS 
BEING  USED. 


ON  OUTPUT: 
G 


REAL(NRD.NN) 

MATRIX  OF  THE  COMPRESSED  PENCIL  AS  DEFINED  ABOVE; 
REAL(NRD.NN) 

MATRIX  OF  THE  COMPRESSED  PENCIL  AS  DEFINED  ABOVE; 
REAL(NRT.NNPM) 

ORTHOGONAL  COMPRESSION  MATRIX  AS  DEFINED  ABOVE; 


INFO  INTEGER 

ERROR  FLAG  WITH  MEANING  AS  FOLLOWS 
=  0  NORMAL  RETURN 

=  NONZERO  IF  SINGULAR  VALUE  DECOMPOSITION  COULD 
NOT  BE  CALCULATED. 

••••• ALGORITHM  NOTES: 

NONE. 


•••••HISTORY: 

THIS  SUBROUTINE  WAS  WRITTEN  BY  W.F.  ARNOLD,  NAVAL'  WEAPONS 
CENTER,  CODE  35104,  CHINA  LAKE,  CA  93555,  AS  PART  OF  THE 
SOFTWARE  PACKAGE  RICPACK,  SEPTEMBER  1983. 


152 


& 


;\ 

i\ 

>5 

v 


SUBROUTINE  RINV(NR,NRD,N,NN,M,E,A,B,CQC,RI 

,G,F,WK1,WRK, RDFLG, 

c 

X 

RFLAG , EFLAG , IBAL , TYPE ) 

Vtf 

c 

•••••PARAMETERS: 
INTEGER  NR, NRD, N 

,NN ,M, IBAL 

CHARACTER  RDFLG, 

RFLAG. EFLAG 

DOUBLE  PRECISION 

E<NR,N) ,A(NR,N) ,B(NR,M) ,CQC(NR,N) , RI(NR,M) , 

X  G(NRD,NN ) ,F(NRD,NN ) ,WK1 (NR,N) ,WRK(N) 

LOGICAL  TYPE 

c 

•••••LOCAL  VARIABLES: 

c 

INTEGER  I.J.K.NPI.NPJ 

c 

c 

n 

NONE. 

c 

•••••SUBROUTINES 

CALLED: 

c 

c 

c 

MULB,  TRNATB 

c 

c 

•••••PURPOSE: 

c 

r 

THIS  SUBROUTINE 

FORMS  THE  MATRIX  PENCIL: 

c 

(  E 

0  )  (A  -B*RI »BT 

> 

c 

LAMBDA*  ( 

)  -  < 

) 

c 

r 

(  0 

ET  )  (  -CQC  -AT 

) 

c 

n 

=:  LAMBDA  «  F 

-  G 

u 

c 

WHERE  T  DENOTES 

MATRIX  TRANSPOSE. 

C  THIS  SUBROUTINE  IS  USEFUL  IN  SOLVING  THE  CONTINUOUS-TIME  GARE. 
C 

C  REF.:  ARNOLD.  W.F.,  "ON  THE  NUMERICAL  SOLUTION  OF 

C  ALGEBRAIC  MATRIX  RICCATI  EQUATIONS,"  PHD  THESIS,  USC, 

C  DECEMBER  1983. 

C 

C  •••••PARAMETER  DESCRIPTION: 

C 

C  ON  INPUT: 

C 

C  NR  INTEGER 

C  ROW  DIMENSION  OF  THE  ARRAYS  CONTAINING  THE  MATRICES 

C  E,  A,  B,  CQC,  RI  AND  WK1  AS  DECLARED  IN  THE  MAIN 

C  CALLING  PROGRAM  DIMENSION  STATEMENT; 

C 
C 
C 
C 


NRD 


INTEGER 

ROW  DIMENSION  OF  THE  ARRAYS  CONTAINING  THE  MATRICES 
F  AND  G  AS  DECLARED  IN  THE  MAIN  CALLING  PROGRAM 
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C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


DIMENSION  STATEMENT; 

N  INTEGER 

ORDER  OF  THE  SQUARE  MATRICES  E,  A  AND  CQC 
ROW  DIMENSION  OF  THE  MATRIX  B; 

NN  INTEGER 

SIZE  OF  THE  MATRIX  PENCIL; 

M  INTEGER 

ORDER  OF  THE  SQUARE  MATRIX  RI 
COLUMN  DIMENSION  OF  THE  MATRIX  B; 

E  REAL (NR, N) 

MODEL  DESCRIPTOR  MATRIX; 

A  REAL (NR, N) 

=  A  -  B*RI*ST  IN  THE  GENERALIZED  CASE; 

B  REAL (NR, M) 

MODEL  INPUT  MATRIX; 

CQC  REAL (NR, N) 

=  CT#Q#C  -  S*RI*ST  IN  THE  GENERALIZED  CASE; 

RI  REAL(NR,M) 

INVERSE  OF  THE  CONTROL  WEIGHTING  MATRIX; 

WK1  REAL (NR, N) 

SCRATCH  ARRAY  OF  SIZE  AT  LEAST  M  BY  N; 

WRK  REAL(N) 

WORKING  VECTOR  OF  SIZE  AT  LEAST  N; 

RDFLG  CHARACTER 

FLAG  SET  TO  ’Y’  IF  RI  IS  A  DIAGONAL  MATRIX; 

RFLAG  CHARACTER 

FLAG  SET  TO  'Y*  IF  RI  IS  OTHER  THAN  THE  IDENTITY 
MATRIX; 

EFLAG  CHARACTER 

FLAG  SET  TO  »Y'  IF  E  IS  OTHER  THAN  THE  .IDENTITY 
MATRIX; 

I3AL  INTEGER 

=  1  IF  CO-ORDINATE  BALANCING  IS  BEING  USED; 

TYPE  LOGICAL 
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£4 


=  .TRUE.  FOR  CONTINUOUS-TIME  SYSTEM 
=  .FALSE.  FOR  DISCRETE -TIME  SYSTEM. 


ON  OUTPUT: 


REAL(NRD.NN) 

PENCIL  MATRIX  AS  DEFINED  ABOVE; 


REAL(NRD.NN) 

PENCIL  MATRIX  AS  DEFINED  ABOVE. 


•••••ALGORITHM  NOTES: 
NONE. 


•••••HISTORY: 

THIS  SUBROUTINE  WAS  WRITTEN  BY  W.F.  ARNOLD,  NAVAL  WEAPONS 
CENTER,  CODE  35104.  CHINA  LAKE,  CA  93555,  AS  PART  OF  THE 
SOFTWARE  PACKAGE  RICPACK,  SEPTEMBER  1983. 
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C 
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C 
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C 

C 

c 

c 
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c 

c 

c 
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c 
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c 

c 
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SUBROUTINE  RESID(NR,NRR,NRW.NRX,N,M,E,A,B,CQC,R,S,RI,RSDM,X,W1, 
X  W2 , WK , IPVT , RTOL , EFLAG , RFLAG , RDFLG , SFLAG , RSD , 

X  TYPE , NOUT ) 


•••••PARAMETERS: 

INTEGER  NR , NRR , NRW , NRX , N , M , IPVT (N ) , NOUT 
CHARACTER  EFLAG , RFLAG , RDFLG , SFLAG 

DOUBLE  PRECISION  E(NR,N) ,A(NR,N) ,B(NR,M) ,CQC(NR,N) ,R(NR,M) , 

X  S(NR.M) ,RI(NR,M) ,RSDM(NRR,N) ,X(NRX,N) ,W1(NRW,N) ,W2(NRW,N) . 
X  WK(N), RTOL, RSD 

LOGICAL  TYPE 


•••••LOCAL  VARIABLES: 
INTEGER  I 

DOUBLE  PRECISION  COND 


•••••FORTRAN  FUNCTIONS: 
NONE. 


•••••FUNCTION  SUBPROGRAMS: 
DOUBLE  PRECISION  D1NRM 


•••••SUBROUTINES  CALLED: 

MADD,  MLINEQ,  MMUL,  MQF,  MSUB,  MULA,  TRNATB 


•••••PURPOSE: 

THIS  SUBROUTINE  CALCULATES  THE  RESIDUAL  MATRIX  AND  ITS  1-NORM 
FOR  THE  GARE  AS  FOLLOWS: 


CONTINUOUS: 

RSDM  =  AT*X*E  +  ET»X»A  -  <BT»X»E+ST)T»RI»( BT»X»E+ST)  +  CQC 


DISCRETE: 


RSDM  =  AT*X*A  -  ET»X»E  ♦  CQC 

-  (BT*X*A+ST)T*( (R+BT*X*B)**-1 )•( BT*X*A+ST) 


WHERE  T  DENOTES  MATRIX  TRANSPOSE. 


REF.:  ARNOLD,  W.F.,  "ON  THE  NUMERICAL  SOLUTION  OF 
ALGEBRAIC  MATRIX  RICCATI  EQUATIONS,"  PHD  THESIS,  USC, 
DECEMBER  1983. 


•••••PARAMETER  DESCRIPTION: 
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c 
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ON  INPUT: 
NR 


NRR 


NRW 


NRX 


CQC 


R 

S 
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INTEGER 

ROW  DIMENSION  OF  THE  ARRAYS  CONTAINING  THE  MATRICES 
E,  A,  B,  CQC,  R,  RI  AND  S  AS  DECLARED  IN  THE  MAIN 
CALLING  PROGRAM  DIMENSION  STATEMENT; 


INTEGER 

ROW  DIMENSION  OF  THE  ARRAY  CONTAINING  THE  MATRIX 
RSDM  AS  DECLARED  IN  THE  MAIN  CALLING  PROGRAM 
DIMENSION  STATEMENT; 


INTEGER 

ROW  DIMENSION  OF  THE  ARRAYS  CONTAINING  THE  MATRICES 
W1  AND  W2  AS  DECLARED  IN  THE  MAIN  CALLING  PROGRAM 
DIMENSION  STATEMENT; 


INTEGER 

ROW  DIMENSION  OF  THE  ARRAY  CONTAINING  THE  MATRIX  X 
AS  DECLARED  IN  THE  MAIN  CALLING  PROGRAM  DIMENSION 
STATEMENT; 


INTEGER 

ORDER  OF  THE  SQUARE  MATRICES  E,  A.  CQC,  RSDM  AND  X 
ROW  DIMENSION  OF  THE  MATRICES  B  AND  S; 


INTEGER 

ORDER  OF  THE  SQUARE  MATRICES  R  »ND  RI 
COLUMN  DIMENSION  OF  THE  MATRICES  B  AND  S; 


REAL (NR, N) 

MODEL  DESCRIPTOR  MATRIX; 


REAL(NR,N) 

MODEL  SYSTEM  MATRIX; 


REAL(NR,M) 

MODEL  INPUT  MATRIX; 


REAL (NR, N) 

MATRIX  PRODUCT  CT»Q*C  WHERE  T  DENOTES  MATRIX 
TRANSPOSE; 


REAL(NR.M) 

CONTROL  WEIGHTING  MATRIX; 


REAL (NR, M) 

STATE  -  INPUT  CROSS -WEIGHTING  MATRIX; 
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RI  REAL (NR, M) 

INVERSE  OF  THE  CONTROL  WEIGHTING  MATRIX; 

X  REAL(NRX.N) 

SOLUTION  MATRIX  FOR  THE  GARE  WHOSE  RESIDUAL  IS  TO  BE 
DETERMINED; 

Wt,W2  REAL(NRW.H) 

SCRATCH  ARRAYS  OF  SIZE  AT  LEAST  M  BY  N; 

WK  REAL(N) 

WORK  VECTOR  OF  LENGTH  AT  LEAST  N; 

IPVT  INTEGER (N) 

WORK  VECTOR  OF  LENGTH  AT  LEAST  N; 

RTOL  REAL 

TOLERANCE  ON  THE  CONDITION  ESTIMATE  OF  R+BT»X*B  WITH 
RESPECT  TO  INVERSION  (DISCRETE  PROBLEM) .  IF  THIS 
TOLERANCE  IS  EXCEEDED  AN  ERROR  MESSAGE  IS  PRINTED 
AND  AN  ERROR  RETURN  IS  MADE; 

EFLAG  CHARACTER 

FLAG  SET  TO  ’Y’  IF  E  IS  OTHER  THAN  THE  IDENTITY 
MATRIX; 

RFLAG  CHARACTER 

FLAG  SET  TO  'Y'  IF  R  IS  OTHER  THAN  THE  IDENTITY 
MATRIX; 

RDF LG  CHARACTER 

FLAG  SET  TO  *Y'  IF  R  IS  A  DIAGONAL  MATRIX; 

SFLAG  CHARACTER 

FLAG  SET  TO  fY'  IF  S  IS  OTHER  THAN  THE  ZERO  MATRIX; 
TYPE  LOGICAL 

=  .TRUE.  FOR  CONTINUOUS-TIME  SYSTEM 
=  .FALSE.  FOR  DISCRETE-TIME  SYSTEM; 

NOUT  INTEGER 

UNIT  NUMBER  OF  OUTPUT  DEVICE  FOR  ERROR  WARNING 
MESSAGES. 

ON  OUTPUT: 


RSDM  REAL(NRR.N) 

THE  RESIDUAL  MATRIX  CALCULATED  AS  INDICATED  ABOVE; 
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IPVT(I)  ERROR  FUG  AS  FOLLOWS 
=  0  NORMAL  RETURN 

s  1  THE  CONDITION  ESTIMATE  OF  R+BT»X«B  WITH  RESPECT 
TO  INVERSION  EXCEEDS  THE  TOLERANCE  RTOL; 


THE  1-NORM  OF  THE  RESIDUAL  MATRIX  DIVIDED  BY  THE 
1-NORM  OF  THE  SOLUTION  MATRIX. 

•••••ALGORITHM  NOTES: 

NONE. 

•••••HISTORY: 

THIS  SUBROUTINE  WAS  WRITTEN  BY  W.F.  ARNOLD,  NAVAL  WEAPONS 
CENTER,  CODE  3510M,  CHINA  UKE,  CA  93555,  AS  PART  OF  THE 
SOFTWARE  PACKAGE  RICPACK,  SEPTEMBER  1983. 
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SUBROUTINE  BALGEN  (N , MA , A , MB , B , LOW , IGH , CSCALE , CPERM , WK ) 

•••••PARAMETERS: 

INTEGER  N,MA,MB,LOW, IGH 

DOUBLE  PRECISION  A(MA,N) ,B(MB,N) ,CSCALE(N) tCPERM(N)  ,WR(N,6) 


•••••LOCAL  VARIABLES: 
NONE. 

•••••FUNCTIONS: 

NONE. 


•••••SUBROUTINES  CALLED: 
REDUCE,  SCALEG,  GRADEQ 


•••••PURPOSE: 

THIS  SUBROUTINE  BALANCES  THE  MATRICES  A  AND  B  TO  IMPROVE  THE 
ACCURACY  OF  COMPUTING  THE  EIGENSYSTEM  OF  THE  GENERALIZED 
EIGEN PROBLEM  A»X  a  (LAMBDA >»B»X.  THE  ALGORITHM  IS  SPECIFI¬ 
CALLY  DESIGNED  TO  PRECEDE  QZ  TYPE  ALGORITHMS,  BUT  IMPROVED 
PERFORMANCE  IS  EXPECTED  FROM  MOST  EIGENSYSTEM  SOLVERS. 

REF.:  WARD,  R.  C.,  BALANCING  THE  GENERALIZED  EIGENVALUE 
PROBLEM,  SIAM  J.  SCI.  STAT.  COMPUT.,  VOL.  2,  NO.  2,  JUNE  1981, 
141-152. 

•••••PARAMETER  DESCRIPTION: 

ON  INPUT: 

MA.MB  INTEGER 

ROW  DIMENSIONS  OF  THE  ARRAYS  CONTAINING  MATRICES 
A  AND  B  RESPECTIVELY,  AS  DECLARED  IN  THE  MAIN 
CALLING  PROGRAM  DIMENSION  STATEMENT; 

N  INTEGER 

ORDER  OF  THE  MATRICES  A  AND  B; 

A  REAL(MA.N) 

CONTAINS  THE  A  MATRIX  OF  THE  GENERALIZED 
EIGENPROBLEM  DEFINED  ABOVE; 

B  REAL (MB, N) 

CONTAINS  THE  B  MATRIX  OF  THE  GENERALIZED 
EIGENPROBLEM  DEFINED  ABOVE; 

WK  REAL (N, 6) 

WORK  ARRAY  THAT  MUST  CONTAIN  AT  LEAST  6»N  STORAGE 
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LOCATIONS.  WK  IS  ALTERED  BY  THIS  SUBROUTINE. 

ON  OUTPUT: 

A.B  CONTAIN  THE  BALANCED  A  AND  B  MATRICES; 

LOW  INTEGER 

BEGINNING  INDEX  OF  THE  SUBMATRICES  OF  A  AND  B 
CONTAINING  THE  NON-ISOLATED  EIGENVALUES; 

IGH  INTEGER 

ENDING  INDEX  OF  THE  SUBMATRICES  OF  A  AND  B 
CONTAINING  THE  NON-ISOLATED  EIGENVALUES.  IF 
IGH  =  1  (LOW  s  1  ALSO),  THE  A  AND  B  MATRICES  HAVE 
BEEN  PERMUTED  INTO  UPPER  TRIANGULAR  FORM  AND  HAVE 
NOT  BEEN  BALANCED; 

CSCALE  REAL(N) 

CONTAINS  THE  EXPONENTS  OF  THE  COLUMN  SCALING  FACTORS 
IN  ITS  LOW  THROUGH  IGH  LOCATIONS  AND  THE  REDUCING 
COLUMN  PERMUTATIONS  IN  ITS  FIRST  LOW-1  AND  ITS 
IGH+1  THROUGH  N  LOCATIONS; 

CPERM  REAL(N) 

CONTAINS  THE  COLUMN  PERMUTATIONS  APPLIED  IN  GRADING 
THE  A  AND  B  SUBMATRICES  IN  ITS  LOW  THROUGH  IGH 
LOCATIONS; 

WK  CONTAINS  THE  EXPONENTS  OF  THE  ROW  SCALING  FACTORS 
IN  ITS  LOW  THROUGH  IGH  LOCATIONS,  THE  REDUCING  ROW 
PERMUTATIONS  IN  ITS  FIRST  LOW-1  AND  ITS  IGH+1 
THROUGH  N  LOCATIONS,  AND  THE  ROW  PERMUTATIONS 
APPLIED  IN  GRADING  THE  A  AND  B  SUBMATRICES  IN  ITS 
N+LOW  THROUGH  N+IGH  LOCATIONS. 

•***• ALGORITHM  NOTES: 

NONE. 

•••••HISTORY: 

WRITTEN  BY  R.  C.  WARD . 
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SUBROUTINE  BALGBK  <N,MZ,Z,M,LOW,IGH,CSCALE,CPERM) 


•••••PARAMETERS: 

INTEGER  N,MZ,M,LOW, IGH 

DOUBLE  PRECISION  Z(MZ,N) ,CSCALE(N) ,CPERM(N) 


•••••LOCAL  VARIABLES: 
NONE. 


•••••FUNCTIONS: 

NONE. 


•••••SUBROUTINES  CALLED: 
GRADBK ,  SCALBK 


•••••PURPOSE: 

THIS  SUBROUTINE  BACK  TRANSFORMS  THE  EIGENVECTORS  OF  A 
GENERALIZED  EIGENVALUE  PROBLEM  A»X  =  ( LAMBDA )*B«X,  THAT  WAS 
BALANCED  BY  SUBROUTINE  BALGEN,  TO  THOSE  OF  THE  ORIGINAL 
PROBLEM. 

REF.:  WARD,  R.  C.f  BALANCING  THE  GENERALIZED  EIGENVALUE 
PROBLEM,  SIAM  J.  SCI.  STAT.  COMPUT.,  VOL.  2,  NO.  2,  JUNE  1981, 
141-152. 


•••••PARAMETER  DESCRIPTION: 
ON  INPUT: 


MZ 


INTEGER 

ROW  DIMENSION  OF  THE  ARRAY  Z  AS  SPECIFIED  IN  THE 
MAIN  CALLING  PROGRAM  DIMENSION  STATEMENT; 


INTEGER 

ORDER  OF  THE  MATRICES  A  AND  B  IN  THE  EIGENPROBLEM ; 


INTEGER 

SPECIFIES  THE  NUMBER  OF  EIGENVECTORS  TO  BE  TRANS¬ 
FORMED; 


REAL(MZ.N) 

CONTAINS  THE  EIGENVECTORS  TO  BE  TRANSFORMED; 


LOW 


INTEGER 

SPECIFIES  THE  BEGINNING  INDEX  OF  THE  SUBMATRICES  OF 
A  AND  B  WHICH  WERE  BALANCED; 


IGH 


INTEGER 
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SPECIFIES  THE  ENDING  INDEX  OF  THE  SUBMATRICES  OF 
A  AND  B  WHICH  WERE  BALANCED; 


CSCALE  REAL(N) 

CONTAINS  THE  REDUCING  COLUMN  PERMUTATIONS  AND 
SCALING  INFORMATION  AS  RETURNED  FROM  BALGEN; 


CPERM 


REAL(N) 

CONTAINS  IN  ITS  LOW  THROUGH  IGH  LOCATIONS  THE 
COLUMN  PERMUTATIONS..APPLIED  IN  GRADING  THE  A 
AND  B  SUBMATRICES  AS  RETURNED  FROM  BALGEN. 


ON  OUTPUT; 

Z  CONTAINS  THE  TRANSFORMED  EIGENVECTORS. 


•••••ALGORITHM  NOTES: 
NONE. 


•••••HISTORY: 

WRITTEN  BY  R.  C.  WARD. 
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